Fish Assemblages of a Sub-Arctic Fjord Show Early Signals of Climate Change Response Contrary to the Benthic Assemblages

Introduction

Climate change has marked effects in the Arctic region: over the last 50 years, it has warmed three times faster than the rest of the world, the extent of sea ice decreased by 43%, and rainfall increased by 24% which increased freshwater runoff and surface stratification (Arctic Monitoring Assessment Program, 2021). Arctic coastal areas, in particular, offer important ecosystem services supporting aquaculture and fisheries, and mitigate global climate change because of high productivity and carbon storage (Johnsen et al., 2009; Smith et al., 2015; Ortega et al., 2019) as well as providing suitable areas for spawning grounds and nurseries for commercially important species (Olsen et al., 2010).

Fjords systems are located at the interface of terrestrial, oceanic, cryospheric, and atmospheric fields and are highly vulnerable to climate change as they have to cope with the cumulative changes that each of these fields is facing (Białogrodzka et al., 2018; Bianchi et al., 2020). A key variable to study climate change impacts from the terrestrial influence is coastal light attenuation. Terrestrially derived organic matters discharge has increased in recent years along the Norwegian coastline due to cross effects between increased precipitations, reduced sulfur deposition, and land-use change (Frigstad et al., 2020). The resulting light attenuation is susceptible to affect photosynthetic organisms, favor tactile predators (e.g., jellyfish) over visual predators (e.g., fish), and finally modify the whole community (Aksnes et al., 2009). On the other hand, changes in the open ocean conditions can also affect the dynamics of species inside the fjords, notably through coastal wind changes resulting in upwelling or downwelling and the vertical stratification at the coast (Svendsen, 1995; Asplin et al., 1999). For instance, an episodic change in wind direction along the Norwegian coast led to the upwelling of cold waters that induced a slowdown in the growth of Pecten maximus inside the fjords (Jolivet et al., 2015). Changes in ocean temperature also cause key predators of sea urchins to fluctuate in population size (Fagerli et al., 2014) which had then consequences on the whole ecosystem due to the grazing pressure sea urchin exert on kelp forest (Strand et al., 2020). Fauna from the open ocean can also enter the fjord and modify its taxonomic composition. For instance, local fjord conditions may attract open ocean megafauna such as humpback whales (Megaptera novaeangliae) and killer whales (Orcinus orca) feeding on wintering Norwegian spring-spawning herring (Clupea harengus) (Jourdain and Vongraven, 2017).

In the south-eastern Barents Sea, 2012 and 2016 are the two warmest years since 1950 regarding surface temperatures while sea ice cover in December 2018 was the lowest since 1951 (González-Pola et al., 2019). These physical changes have already resulted in biodiversity changes in the Barents Sea and induced a poleward shift of fish assemblages and an increase of warm water boreal species (Wassmann et al., 2011; Fossheim et al., 2015; Frainer et al., 2017). Similarly, for megabenthic species, a recent poleward expansion of boreal assemblages and a corresponding decrease in the importance of arctic assemblages were reported (Jørgensen et al., 2019). These changes in megabenthic assemblages are also visible in arctic fjords such as Kongsfjord in Svalbard (Kedra et al., 2010b) where it has eventually led to a regime shift of the rocky-bottom assemblages with a sharp increase of macroalgae cover (Kortsch et al., 2012). Megabenthic species are considered as a good indicator of long-term changes in climatic conditions since they are largely non-mobile and long-lived, and thus representative of the sampling site (Reynoldson and Metcalfe-Smith, 1992; Kedra et al., 2010b; Blicher et al., 2015), while fish species, more mobile, are hypothesized to react more quickly to climatic changes (Kortsch et al., 2015). It is thus interesting to compare fish and benthic assemblage trends facing climate change as fish may provide early warning of the changes that might later also affect benthic assemblages.

Climate change is not the only threat to fjord ecosystems. In link with increased anthropogenic activities, 11 non-indigenous species were introduced in the Barents Sea region over 1960–2015, among which 55% got established (Chan et al., 2019). In fjords of Northern Norway, fish and megabenthos assemblages have to face invasion by red king crab Paralithodes camtschaticus, which is one of the rare large, higher-trophic-level organisms that has established itself into new areas (Jamieson et al., 1998; Falk-Petersen et al., 2011). Red king crabs were introduced from the Pacific to the Russian Barents Sea at the end of the 1960s, and reached the Norwegian coast in the early 1990s (Oug et al., 2018). Those are active predators (Jørgensen and Primicerio, 2007) that led to substantial structural and functional changes in megabenthic assemblages with the decrease of large suspension and surface deposit feeders in adjacent fjords while small, mobile, shallow burrowing and predator species with planktonic larval development increased (Oug et al., 2018). Red king crabs also prey on benthic fish eggs (e.g., lumpfish and sculpin) but are of minor importance in the diet of large demersal fish and cod that prey on small and medium red king crabs (Pedersen et al., 2018).

Fjords can harbor relict population of cold-water species, as it has been shown with the presence of the cold-water species Bathyarca glacialis in the Olsofjord during the warming that occurred in the younger Dryas (12,000 years ago) (Spjeldnæs, 1978). Today, sub-arctic fjords of Northern Norway are located at a latitude where warm-water boreal taxa dominate offshore (Johannesen et al., 2012; Jørgensen et al., 2015). Very few fjords have an inner basin isolated from the rest of the fjord and open ocean by a sill, in which water remains very cold (in Norway, it is essentially Porsangerfjord, Balsfjord, and Ullsfjord; T. Pedersen (UiT) personal communication). Such fjords harbor a local and probably isolated population of cold-water arctic species, as it has been observed in Porsanger (Oug and Fuhrmann, 2013; Fuhrmann et al., 2015), and could act as a refugee area for these species (Włodarska-Kowalczuk et al., 1998; Drewnik et al., 2017; Yletyinen, 2019). Those fjords can be considered as a “mini arctic” (Oug and Fuhrmann, 2013), providing unique insights into the changes undergone by arctic assemblages facing climate change and invasive species, that might be of relevance for Svalbard and Northern Barents Sea. Here we focus on the Porsangerfjord as a case study of the dynamics of the biodiversity of sub-arctic fjords facing multiple stressors, i.e., climate change and invasive species. This fjord was closed to shrimp fishing in the early 1970s after intensive fishing caused the overexploitation of cod and small young fishes (Søvik et al., 2020). Hence, it is considered to be free of bottom fishing impacts. The fjord also harbors a large population of red king crab, that entered the fjord in the mid-1990s, and progressed toward the inner-most part of the fjord that they reached after 2010 (Fuhrmann et al., 2015; Pedersen et al., 2018). It is now an important commercial red king crab fishery with 921 tones caught in the Porsangerfjord in 2018 (Søvik et al., 2020). Climate change, through increased water temperature and light attenuation, has been recorded in the fjord. Over 1986–2015, air temperature increased by 0.05°C/year near the fjord mouth and 0.04°C/year in the inner basin (Cieszyńska and Stramska, 2018). Over 2009–2014, light attenuation is estimated to be quite high in the inner part of the Porsanger (up to K_d(PAR) = 1.96, Trine Bekkby (NIVA) unpublished data, Frigstad, 2020). Light attenuation as a function of increased precipitation was approximated here by including riverine discharge among the environmental variables used in the study. More generally, depth, sediment type, temperature, current, and sea ice have been related to variation in the structure of megabenthic assemblages in the Arctic and their effects were thus investigated here (Gray, 2002; Piepenburg, 2005; Sejr et al., 2009; O’Sadnick et al., 2020). For fish assemblages, depth and temperature gradients appear to be the main structuring variables in the Barents Sea (Johannesen et al., 2012), while sediment type and seabed currents are also reported to be important structuring variables of bentho-demersal fish assemblages (Vaz et al., 2007; Pickens et al., 2021). Since there is a strong salinity gradient in the Porsanger, salinity was also included as a potential driver of assemblages in this study for both fish and megabenthos.

Here we aim at investigating the consequences of climate change and the increase in red king crab biomass on fish and megabenthic assemblages. More specifically, we investigated which of these pressures are the most structuring variables for fish and megabenthic assemblages, whether their importance changed over 2007–2019, and if the mobile fish assemblages reacted earlier than the less mobile megabenthic assemblages. Two datasets covering 2007–2019 for fish and megabenthic species respectively were explored as functions of environmental variables known to influence fish and megabenthic species distribution in fjords. We first investigated (1) temporal trends in the biomass of each fish and megabenthic species and (2) temporal changes in biodiversity indices. Then, we characterized (3) whether a shift is visible in the relationship between species and environmental variables (i.e., a dramatic change in species composition or environment); (4) the average structure of the association of fish and megabenthos with environmental variables (the distinction between cold and warm water species) and (5) the interannual variations or the importance of each environmental variable.

Materials and Methods

Study Area

The Porsanger, located in Finnmark above the polar circle is one of Norway’s largest fjord. It is 123 km long and covers an area of about 1,877 km². This large and open sub-arctic fjord is characterized by the heterogeneity of its environmental conditions with a strong along-fjord gradient (Fuhrmann et al., 2015). Three areas can be distinguished in the fjord: (1) the outer fjord faces the Barents Sea and is influenced by warm and saline coastal waters with depths around 300 m; (2) the center of the fjord is separated from the outer fjord by a sill that is about 200 m deep: (3) the inner fjord is isolated from the center of the fjord by a narrowing of the fjord width and a sill that is about 60 m deep. In the inner fjord, seabed temperatures are below zero for half of the year (Fuhrmann et al., 2015) and freshwater inputs from the three rivers in the inner fjord make the salinity lower than in the rest of the fjord (Myksvoll et al., 2012). This leads to a sea ice coverage, present every year but varying in extent, during late winter until the beginning of the spring (Megan O’Sadnick, SINTEF, pers. comm.). This creates special conditions that enable the presence of unique species compared to the adjacent fjords at the same latitude (Søvik et al., 2020). Freshwater discharge and stratification in the Porsanger are low (Svendsen, 1995). There is an anticlockwise surface water circulation in the fjord with high saline and warm surface waters entering the fjord on the west and low saline water exiting the fjord on the east (Svendsen, 1991).

Environmental Variables

Depth data were downloaded from the GEBCO website¹ and sediments were taken from the ‘‘Marine---Seabed sediments (grain size) detailed’’ product from the Geological Survey of Norway (NGU) website.² The 17 categories of sediments originally documented in the sediment data were re-coded into 3 main categories to avoid loss of power in the analysis: mud, sand, and rough substratum. Surface and bottom temperatures, salinity as well as seabed current were extracted from the NordKyst800 model outputs at 800 × 800 m resolution (Albretsen et al., 2011; Asplin et al., 2020). An evaluation of this model was performed (Supplementary Figure 1). The mean errors between observations and model do not have significant temporal trends, but they are significantly different between areas. Despite this spatial bias, the along-fjord gradient was well reconstructed by the model. Outputs of the NorKyst800 models could thus be used to identify potential changes in the relationship between environment and taxa along the fjord. River discharges (m³/s) were taken from the outputs of an HBV (Hydrologiska Byråns Vattenbalansavdelning) model run by the Norwegian Water Resources and Energy Directorate (NVE) (Bergström, 1995; Lindström et al., 1997). We selected the cells corresponding to the inner and central/outer fjord and compute the annual mean. Sea ice maximum extent was obtained from automatic analysis of Moderate Resolution Imaging Spectroradiometer (MODIS) images (O’Sadnick et al., 2020). Sea ice is only present in the inner part of the fjord, so NAN was attributed to sampling sites in the central and outer fjord. Due to the annual or less than the annual frequency of our samplings of a large number of taxa with different life cycles covered in this study, we used annual means for the environmental variables that are used to potentially explain shifts in the assemblage composition.

To identify potential temporal trends in surface temperature, bottom temperature, surface salinity, bottom salinity, and seabed currents a linear regression of the values over time was applied in each pixel (800 × 800 m) (Figure 1). The slope of the regression was printed only in the pixels with a p-value < 0.1. For non-spatialized variables (sea ice maximum extent and river runoff) similar linear regressions were computed independent of pixels. For visualization, a GAM model was also fitted.

Figure 1. Temporal trends of the environmental variables. For spatialized data (A–E), linear regression was computed in each 800 × 800 m cell, only slopes with p-values < 0.1 are printed. For non-spatialized data (F–H), linear regression was performed on the annual mean (stratified biomass), for the inner fjord for sea ice extent (F) and separated between inner (G) and central/outer (H) areas for river runoff.

Biological Sampling

Benthic assemblages were sampled in 2007, 2008, 2009, 2010, 2014, 2016, and 2019. Each of these years, between 5 and 22 sites in total were visited at the outer, center, and inner parts of the fjord (Figure 2). Hauls were conducted for 5 min at a mean speed of 1.5 knots with a 2 m beam trawl (codend 4 mm) and washed through a 5 mm sieve. The catch was then sorted, and each taxon was counted and weighted (i.e., biomass is wet mass). Megabenthic taxa were identified to the lowest possible taxonomic level. However, as this was not always possible, some taxa were nested (a species and its genus were present in the database) (see Supplementary Table 1). Fish assemblages in the Porsanger were sampled as part of the annual coastal ecosystem survey in October. During 2007–2019, this survey sampled every year between 2 and 8 sites in the fjord (Figure 2). Due to this low number of samplings and since fish species are mobile in the fjord, outer and central area were pooled for fish data. The year 2008, which has only 2 sampling sites into the Porsanger was removed from the dataset. Fish and megabenthic species (> 1 cm, Stratmann et al., 2020) were collected with a Campelen 1,800 bottom trawl, rigged with rockhopper ground gear, and towed on double warps. Sampling effort was standardized to a fixed effort of 15 min [equivalent to a towing distance of 0.75 nautical miles (∼1.4 km)]. The horizontal opening of this trawl is 17 m, and the vertical opening is 4 m (Engås and Ona, 1990). The mesh size is 80 mm (stretched) in the front and 24 mm at the codend, allowing the capture and retention of small fish and the largest megabenthos (benthic megafauna larger than 4 mm) from the seabed. This represents a total of 78 valid hauls.

Figure 2. Sampling sites for megabenthic taxa (A) and fish taxa (B) in the inner (green), central (yellow), and outer (red) parts of the fjord.

Since the number of sampling sites varied between years (Supplementary Figures 2A,C), a non-parametric Kruskal-Wallis test was performed on each of the megabenthos and fish datasets to determine if there was a significant difference in the number of sites sampled between years. As the number of sampling sites is low in some years, the power of the test was also calculated with the kwpower function from the MultNonParam R package (Kolassa and Jankowski, 2021).

Some species are not always correctly identified and were grouped to a higher taxonomic level (i.e., Benthosema glacialis, Myctophidae, and Myctophiformes under Myctophiformes; Argentina silus and Argentina shyraena under Argentina; Sebastes mentella, Sebastes norvegicus, Sebastes viviparous, and Sebastes under Sebastes). For both fish and megabenthic taxa, to avoid bias caused by sporadic catches of taxa that can be due to a restricted spatial and temporal coverage of our data and/or the non-homogenous accuracy in taxonomic identification across different crews, only taxa occurring in more than 5% of all sampling hauls and over more than one-third of the time series were kept for analysis. A total of 24 fish and 71 megabenthic taxa were selected (the full list is given in Supplementary Table 1) and their biomass was standardized by the area sampled during each haul (computed from the distance trawled and the width of the trawl). Biomass was log(x + 1)-transformed to reduce variability caused by the catch of numerous individuals of the same taxon that may occur for some of the fish schools or sessile megabenthic taxa.

Temporal Trends in Biomass of the Taxa and Biodiversity Indices

The trend of the biomass of each fish and megabenthic taxa was computed as a function of time with a regression, in each area of the fjord (inner, central, outer for megabenthic taxa and whole fjords for fish taxa). Since the time series are not very long and we are interested in trends, a simple linear regression was computed. The stratified biomass in each area was computed to account for the different number of sampling sites in each area and each year (i.e., mean biomasses for each taxa in each area and each year were multiplied by the surface of the area and divided by the total surface of the fjord, these weighed mean biomass were then summed to get the total biomass of the whole fjord, Bellail et al., 2001). For megabenthos, since each area was considered separately, the stratified biomass was the mean biomass in each area each year. Biomass was log(x + 1)-transformed to reduce variability. Finally, the regression was performed on strictly positive biomass only, as it is impossible to determine whether a taxon was really absent in an area a specific year or if we failed at sampling it.

Diversity indices were computed for each area for megabenthos and for the whole fjord for fish data to characterize the dynamics of taxa composition and the changes undergone in the fjord in recent years. In each area for the megabenthos and the whole fjord for megabenthos and fish, log(x + 1)-transformed stratified biomass, log(x + 1)-transformed stratified abundance, taxonomic richness, Shannon-Weiner diversity index that informs on community diversity and Pielou’s evenness that informs on the equality of biomass distribution between taxa (Pielou, 1966) were computed. Pielou’s evenness is independent of the number of taxa and is only influenced by the distribution of individuals among taxa. Specifically, to follow the evolution of the particular ecosystem in the inner fjord, the proportion of the stratified biomass accounted by each taxon each year was represented and the 5 taxa with the largest stratified biomass were identified each year (see Supplementary Figure 3).

Statistical Analysis: The STATICO Method

The STATICO analysis belongs to the ordination analysis family that enables the study of the relationship between taxa and their environment (Thioulouse et al., 2004). This analysis has the particularity to handle complex relationships between numerous taxa and environmental variables with spatially and temporally resolved information. STATICO is based on a series of k (number of years) times two tables: a site × taxa table, which provides biomass for each taxon at each site sampled, and a sites × environment table, which provides values for each environmental variable at each site. To avoid placing undue emphasis on the taxa with low frequency, Hellinger’s transformation was applied on biomass according to Rao (1995). Some sediment types were sampled many times while others were not. To remove this unbalanced sampling among sediment types (Supplementary Figures 2B,D), sediments were grouped in three categories and coded as a numeric variable with larger values representing larger particle sizes (mud = 1, sand = 2, rough substrate = 3).

The first step of STATICO consists of co-inertia analysis that combines, for each year, the sites × taxa table with the sites × environments table. This yields a cross-covariance table between taxa and environment for each year. A preliminary step to the co-inertia analysis is to perform two normed principal component analyses (within-PCA), to remove the effect of time: one for the sites × taxa table and one for the sites × environments table. During this step, in each of the two tables, variables are standardized separately within each year (Thioulouse et al., 2018). Since this step assumes a linear relationship between taxa’s biomass and environmental variables, the relation between each pair of taxon and environmental variables was visually verified to be roughly linear. After the co-inertia analysis, a Monte-Carlo permutation test of rows and columns of the taxa and environment tables was performed, with the statico.krandtest function from the ade4 R package (Thioulouse et al., 2018), for each year, to test the significance of the co-variations between taxa and environment. The second step of STATICO consists of partial triadic analysis (PTA: Tucker, 1966; Thioulouse et al., 2018) to analyze the series of cross-covariance tables across years.

1. The interstructure identifies similarities in the relationships between taxa biomass and the environment between years and enables to detect potential shift in the time series. For each year with every other year, vectorial correlation coefficients Rv are calculated between cross-covariance (environments × taxa) tables and gathered in a year × year Rv matrix to identify similarities between years. On the correlation circle, the length of one arrow on the first axis represents the weight of the cross-covariance (environments × taxa) table of this year in the computation of the compromise (the averaged cross-covariance table). The larger this length, the greater the contribution of a year to the compromise.

2. The compromise assesses temporally stable relations between taxa biomass and environmental variables. The first eigenvector of the eigen analysis of the Rv matrix is used to weight the cross-covariance tables and calculate their weighted mean. A high weight for a given year indicates that this year is representative of the relation between taxa and the environment during the time series while a low weight indicates an unusual relation this year. Similar weights for all years indicate stable relations between taxa and the environment. The compromise is then calculated using a centered PCA.

3. The intrastructure assesses temporal changes in relations between taxa biomass and environmental variables. To assess potential spatial patterns in taxa–environment relations, we projected the rows (sampling sites) of the initial sites × taxa and sites × environments tables onto the compromise space as supplementary elements (Thioulouse et al., 2018). These projections provide an overview of annual variations in taxa–environment relations, also called trajectories.

One STATICO analysis was run for megabenthic taxa and one for fish taxa. At the compromise step, a K-means clustering algorithm was applied to the taxa in the compromise space to determine groups of taxa with a common relationship with the environment. The number of groups was given by the Calinsky-Harabasz criterion (Calinski and Harabasz, 1974). Taxa with the largest biomass in each of the fjord areas were computed based on the stratified biomass of each taxon each year. Finally, stratified biomasses of taxa belonging, respectively, to each K-means cluster (warm-water and cold-water taxa) were summed and log(x + 1) transformed. A linear regression of these log(x + 1) transformed stratified biomass of each group was fitted to evaluate the temporal trend of their biomass. Finally, at the intrastructure step, to disentangle the cause of the variations in the stable variables (i.e., depth and sediment), a Kruskal-Wallis and a Pearson’s Chi-square test were computed. This enables to disentangle if there was a significant difference in sampling across depth or sediments or if taxa were moving.

Results

Environmental Variables

Over 2007–2019, a significant decrease in surface and bottom salinity is visible in outer, central, and inner parts of the fjord with a loss of 0.01–0.05 psu/year (Figures 1C,D). Seabed current shows some significant increase of 0–0.001 m/s near the coast of the central fjord (Figure 1E). There are no significant trends for surface and bottom temperature (Figures 1A,B). Sea ice maximum extent shows some decreasing trends between 2010 and 2015, but no significant decrease over 2007–2019 (Figure 1F). River runoff into the inner and medium/outer parts of the fjord has increasing, but not significantly (Figures 1G,H).

Temporal Trend of the Biomass of Each Taxon

For megabenthic taxa, the biomass of Norwegian shrimp Pontophilus norvegicus shows a significant increase in the outer fjord while the biomass of brittlestars Ophiuridae significantly decrease in the inner fjord similarly to Island scallop Chlamys and sea anemones Actiniaria in the central fjord (Table 1). In addition, some non-significant trends also appear important to be reported: red king crab Paralithodes camtschaticus have an increase of its biomass in the inner fjord. Regarding fish, the biomass of haddock Melanogramus aeglefinus and Atlantic hookhead sculpin Artediellus atlanticus significantly decrease while the biomass of redfishes Sebastes significantly increase. A non-significant trend is also worth noting, with the decrease in the biomass of Atlantic wolffish Anarhichas lupus in the fjord over 2007–2019 (Table 2).

Table 1. Coefficient and p-value of the linear regression of the log-transformed stratified biomass of megabenthic and fish taxa.

Table 2. Coefficient and p-value of the linear regression of the biodiversity indices (taxonomic richness, log(x + 1) transformed stratified biomass, log(x + 1) transformed stratified abundance, Shannon diversity, and Pielou’s evenness) for megabenthic and fish taxa over time.

Temporal Trend of the Diversity Indices at the Scale of the Assemblage

Taxonomic richness and Shannon diversity index significantly increased over 2007–2019 for the fish dataset (Table 2). For megabenthic taxa, there is a significant decrease of Pielou’s evenness in the whole fjord, which indicates that the biomass is becoming less evenly distributed between taxa, with the emergence of a dominant taxon (Table 2). More specifically, toward the end of the time series (2014 and 2019), the biomass of the red king crab become among the 5 largest biomasses in the inner fjord (Supplementary Figure 3). In the inner fjord, a decrease in the abundance of megabenthos is also visible, but non-significant. Finally, a non-significant increase in the biomass of megabenthic taxa in the central fjord can also be noted (Table 2).

STATICO Analysis

For both megabenthic and fish taxa, the Kruskal-Wallis tests indicated no significant differences in the number of sampling sites between years (for megabenthos: Kruskal-Wallis Chi-squared = 7, p = 0.429, power = 0.999; For fish: Kruskal-Wallis Chi-squared = 12, p = 0.446, power = 0.999).

The co-variation between taxa and environment was significant all years, except in 2010 and 2019 for megabenthic taxa, and in 2009 for fish. This highlights that those years the assemblages were likely driven by other variables than those considered in this study.

The interstructure step of the STATICO analysis quantifies how much the relationship between taxa and environment is similar between years. The first axis of the interstructure accounts for the largest share of the similarity in the taxa-environment relationship between years. For megabenthic taxa, the first axis shows moderate common structure between years (30.0% of explained variability) (Figure 3B), while for fish the first axis explained more of the correlation between years (46.62%) (Figure 3E). In addition, for megabenthic taxa, values of the vectorial correlation coefficient Rv are moderate to low which indicate that relation between taxa and the environment tend to vary between years (Figure 3C). Conversely, for fish, the relation between taxa and the environment was more stable as shown by the moderate to high values of Rv (Figure 3F). Similarly, the length of the arrows on the first axis of the fish correlation circle is long and quite even (Figure 3D), indicating relative stability of the taxa-environment relation through the years. For megabenthic taxa, the length of the arrows on the first axis is less even, pointing again toward the variability of the taxa-environment relation for those taxa (Figure 3A).

Figure 3. Interstructure analysis with correlation circle with scores for years for megabenthos (A) and fish (D). (B,E) Eigenvalues of the Rv matrix that gathers vectorial correlation in taxa-environment covariances between 2 years). (C,F) Histogram of the Rv values between 2 years.

The compromise analysis informs on the average relation between taxa and environment through the time series. For megabenthic taxa, the first axis accounts for 79.1% of the variability in the compromise space, and the second axis 13.5%. For fish, the first axis explains 84.6% of the variability in the compromise space while the second axis accounts for 14.2%. In both cases, the first axis represents the strong structuring importance of the along-fjord gradient, with the outer fjord deep sampling sites with high values of temperature, salinity, and currents on the left and the inner fjord shallow sampling sites with more river runoff, and sea ice on the right (Figures 4B,D). The second axis, in both cases, is associated with sediments, with a gradient of grain size from the bottom to the top of the second axis (mud, then sand, then coarse sediments). For megabenthos, bottom salinity, followed by surface temperature, river runoff, sea ice extent, and bottom temperature and depth are the main drivers, while for fish assemblages it is bottom salinity, then depth, surface salinity and surface temperature.

Figure 4. Compromise analysis for megabenthos and fish. (A,C) The space of the taxa and (B,D) the space of the environmental variables.

The clustering of the taxa along the scores on the first axis of the compromise lead to two groups for both megabenthos, and fish taxa (Figures 4A,C). Megabenthic taxa can be distinguished between cold water taxa (Figure 4A, cluster 1) associated with high values of sea ice extent and river runoff and warm waters taxa (Figure 4A, cluster 2) associated with high values of temperature, salinity, and current. The largest biomass of cold-water megabenthic taxa are formed by ascidians Ascidiacea, cookies cutter sea star Ctenodiscus crispatus, red king crab Paralithodes camtschaticus, whelks Buccinidae, and soft coral Nephtheidae while the largest biomass of warm water megabenthic taxa originate from shrimps Pandalidae, sea anemones Actiniaria, brittlestars Ophiuridae, Nephtyidae and Sevenline shrimp Sabinea septemcarinata. Similarly, fish taxa can also be partitioned between warm water taxa (Figure 4C, cluster 1) associated with warm and saline deep water with strong current and cold-water taxa (Figure 4D, cluster 2) found in shallower waters with lower values of temperature, salinity, and current.

The biomass of the warm-water and the cold-water megabentic groups did not show any significant temporal trends (Supplementary Figure 5A). It can be noted that when removing red king crab from the cold-water taxa, no significant trends were observed either (Supplementary Figure 5B). The biomass of the cold-water fish taxa significantly decreased in the whole fjord while the biomass of warm-water fish taxa did not vary significantly over time (Figure 5). In addition, the importance of bottom temperature as a structuring variable of the assemblages sharply increased after 2016 (Figure 6B). This variable was the one showing the largest variation in its importance in structuring the assemblages (large size of the ellipse) (Figure 6A). Conversely, for megabenthos the structuring effect of environmental variables appears relatively stable through the period considered (2007–2019) as indicated by the medium size of the ellipses (Supplementary Figures 6A,B). Fish and megabenthos taxa maintained stable relations with the environment through the years (Supplementary Figures 6C–F). For variables that do not change in time, depth and sediments, their variations could be interpreted either because the significantly different depth and sediments were sampled through time or because taxa moved. For depth, there is no significant difference in the sampled values between years (For megabenthos: Kruskal-Wallis Chi-squared = 10.64, p-value = 0.1; for fish: Kruskal-Wallis Chi-squared = 2.39, p-value = 0.967) as well as for sediments (for megabenthos: Chi-squared = 12.653, p = 0.554; for fish: Chi-squared = 2.81, p-value = 0.993) meaning that the observed variability can be caused by taxa moving among different sediment types and depth.

Figure 5. Linear regression of the log(x + 1) transformed stratified biomass of fish warm-water and cold-water taxa (from the K-means clustering on the score of the two first axis of the compromise) over time. Equation of the regression, its r², and its p-value are printed for each group for the whole fjord.

Figure 6. Intrastructure analysis for fish: (A) The yearly average position of environmental variables in the space of the compromise is the center of the ellipses and each point is the yearly position (two last numbers of the year are printed), (B) score of the environmental variables on axis 1 of the compromise.

For each part of the fjord, the 5 megabenthic taxa with the largest stratified biomass were reported: (1) for inner fjord: cookie-cutter seastar Ctenodiscus crispatus, sea urchins Strongylocentrotus, shrimps Pandalidae, sea urchins Echinus, and basket star Gorgonocephalus; (2) for central fjord: red king crab Paralithodes camtschaticus, soft coral Nephtheidae, shrimps Pandalidae, ascidians Ascidiacea and polychaetes worms Maldanidae; (3) for outer fjord: sea anemone Actiniaria, shrimps Pandalidae, hairy cockle Ciliatocardium ciliatum, brittle star Ophiuridae and polychaetes worms Ampharetidae. Similarly, fish taxa with the largest biomass in the whole fjord were haddock Melanogrammus aeglefinus, cod Gadus morhua, American plaice Hippoglossoides platessoides, whitch flounder Glyptocephalus cynoglossus, starry ray Amblyraja radiata.

Discussion

We used a complex multivariate analysis to identify the trends of the co-variation with the environment of megabenthic and fish taxa in a sub-arctic fjord facing both climate change effects and invasion by red king crabs. Our results show that the fish assemblages are increasingly driven by bottom temperature, especially after 2016. Those assemblages are markedly changing at the detriment of the cold-water boreal fish taxa whose biomass significantly decreased over 2007–2019, while the biomass of warm-water fish taxa did not show a significant trend. In particular, we observed a significant decreasing trend for the cold-water haddock Melanogrammus aeglefinus and Atlantic hookhead Artediellus atlanticus. Increased taxonomic richness and Shannon diversity index indicate that warm-water fish species that were occasionally present in the fjord at the beginning of the time series tend to become present every year, toward the end of the time series. This included argentine Argentina, blue whiting Micromesistius poutassou, lanternfishes Myctophiformes, and round ray Rajella fyllae (Supplementary Figure 7).

Fish and megabenthic taxa co-variated along the fjord, discriminating between warm-water taxa present in the outer part of the fjord, and cold-water taxa present in the central and inner part. Such co-variation of fish and megabenthos with environmental gradients has been mapped for the Barents Sea (Johannesen et al., 2017) but without investigating temporal responses to climate change between fish and megabenthic assemblages, as done in our work. Over the studied period, the main change in environmental variables that the taxa had to cope with was the freshening of the water in the fjord and the decrease of the maximum sea ice extent. This density reduction could strengthen water stratification in the fjord, which has led to dissolved oxygen decline in fjords in southwest Norway (Aksnes et al., 2019). However, the rate of renewal of the deep water below the sill in the inner basin (Askeland, personal communication) would likely prevent the Porsanger from future anoxia. This seems confirmed by the O₂ measurements realized in the whole fjord in May 2019, with the lowest values > 6.5 mmol.m^–3 (IMR coastal survey, unpublished data). A decrease in sea ice maximum extent is particularly visible in 2010–2015. Yet, over 2007–2019, no significant decrease in maximum sea ice extent, increase in surface or bottom temperature could be detected in the fjord. However, the warming of waters along the North Norwegian coast due to a larger penetration of warm Atlantic waters (Asbjørnsen et al., 2020; Skagseth et al., 2020) has most likely driven the arrival of warm-water fish species similar to what has been reported in the Barents Sea (Fossheim et al., 2015).

With climate change and increased sea temperature, megabenthic species composition is also hypothesized to change toward more warm-water species (Kedra et al., 2010b). These changes are thought to be particularly strong for arctic benthic species that are documented to have a narrow environmental tolerance (i.e., stenothermal species) and lower upper thermal limit compared to arcto-boreal and boreal species (Morley et al., 2019; Renaud et al., 2019). In the Porsanger, including in the Arctic inner part of the fjord, no significant decrease of the biomass of cold-water species was recorded. However, given that we did not observe a significant increase in bottom temperature inside the fjord, this finding might not be surprising. This lack of a significant trend in megabenthic diversity is in line with results from previous studies. For instance, Węsławski et al. (2017) reported that marked changes in the physical environment had little effect on the benthic fauna of two fjords of the West Spitsbergen, because of the fauna’s adaptation to the scale of the ongoing changes. Indeed, megabenthic assemblages might show some degree of stability facing climate change due to their inherent adaptive capacities to survive long periods of seasonally low food availability (Sun et al., 2009) as well as opportunistic feeding strategy (Iken et al., 2010; Weslawski et al., 2011). Renaud et al. (2019), which modeled megabenthic species habitat under an end-of-century scenario of warming and acidification, found that mean habitat loss was small (0–11%) and that there was no significant difference in habitat loss between arctic and arcto-boreal species. In addition, the response of megabenthic cold-water taxa to climate change sometimes cannot be fully evaluated considering the environment the taxa live in, and notably their thermic preferences. Other parameters not considered in this study, such as metabolic characteristics, that enable species to feed on zooplankton (Sun et al., 2009) or macroalgae detritus (Renaud et al., 2015) instead of on sea-ice dependent micro-algae, might also play a major role in climate change adaptation. Lowering in size could also favor the survival of species as it reduces their metabolic energy demand (Garilli et al., 2015). Further, macroalgae coverage experienced a very sharp decrease due to the burst in sea urchin abundance that happened into the fjord from the beginning of the 1970s following the overexploitation of Norwegian spring-spawning herring, large cod, and saith (Larsen, 2010). Macroalgae offer important habitats for numerous benthic species and juvenile fish (Larsen, 2010). The recovering in the extent of macroalgae might have been a favorable parameter for megabenthic species that might compensate to some extent the changes in environmental variables (decreased salinity, increase bottom current, and decreased sea ice extent). Finally, this lack of significant trends could also arise from bias in our data. The patchiness of the sediments, especially in the south-central fjord, makes it difficult to get good coverage of the different sediment types where to follow a trend analysis in time. Homogenous areas for long-term monitoring should be identified and sampled every year. Further, only species occurring in more than 5% of the sampling sites and one third of the years of the time series were kept in this analysis. Rare species could provide early signals of environmental shifts and their absence from this dataset could also explain the lack of significant trends observed.

Climate change was also hypothesized to lead to more homogenous benthic assemblages (Kedra et al., 2010b; Weslawski et al., 2011). But such homogenization was not observed in the Porsanger since there were no significant decreasing trends of Shannon diversity for megabenthic taxa. In the whole fjord, the significant decrease of Pielou’s evenness for megabenthic taxa highlights that the biomass of some few taxa become dominant through 2007–2019. The increase of red king crab biomass in the inner fjord (non-significant), as well as its biomass becoming amongst the 5 largest in the inner fjord from 2014, points toward red king crab becoming a dominant taxon in the inner fjord. Red king crabs were reported to negatively impact the productivity and the biomass of megabenthic assemblages (Pedersen et al., 2018). Here we observed a decrease in the abundance of megabenthic organisms in the inner fjord, which is a sign of red king crab impact on benthic assemblages that was also reported by Falk-Petersen et al. (2011) in invaded areas in the Barents Sea. The significant decrease of brittle stars Ophuroidae in the inner fjord, as well as the one of Island scallops Chlamys in the medium fjord, could thus originate from the predation of the red king crabs on these taxa (Jørgensen and Primicerio, 2007). Since its introduction, benthic taxa composition has changed in three fjord areas in Northern Norway including Porsanger, with a decrease of echinoderms and most bivalves while small polychaetes prospered (Oug et al., 2018). Some further decrease of in the abundance of benthic taxa, especially in the cold-water assemblages in the inner fjord, are thus expected as the red king crab continues its progression. Red king crabs are estimated to consume between 1 and 18% of the production of benthic invertebrates in the Porsanger (Fuhrmann et al., 2015). Søvik et al. (2020) indeed noted the attractiveness for red king crabs of the inner fjord since those areas were later invaded and therefore probably contain more food. Among megabenthic taxa with the largest biomass in the inner part of the fjord, cookie-cutter sea star Ctenodiscus crispatus and sea urchins Strongylocentrotus are widely distributed through the Barents Sea (Jørgensen et al., 2015). The inner fjord megabenthic assemblage thus appears as a mix between sub-arctic and a few arctic taxa, such as basket stars Gorgonocephalus, Iceland scallop Chlamys islandica, hairy cockle Clinocardium ciliatum, Northern yoldia Yoldia hyperborean, or crangonid shrimp Sclerocrangon ferox. This is in line with the dominance of sub-arctic (arctic-boreal) taxa observed in the arctic fjord Hornsund (Svalbard) (Kedra et al., 2010a). In addition, C. crispatus is a detritivore sea star feeding on small particles that settle on the seabed. This might indicate a weak current regime of the ecosystem with a large quantity of detrital matter that represents C. crispatus food supply and allow it to sustain its large population (Jørgensen et al., 2015). C. crispatus also appears to be a prey of red king crabs (Jørgensen and Primicerio, 2007) with the potential to eradicate it from the area, as it happened in an adjacent fjord (Varangerfjorden) over 1994–2007 (Oug et al., 2011). A major regime shift is thus a potential threat for the inner cold parts of the fjord if red king crabs continue their progression. In addition, a decrease in the condition of red king crabs has been reported by fishermen, at the beginning of 2021, with a large share of the catch composed of small crabs with poor meat content where catches were previously dominated by large crabs with good meat content (H.K. Strand, IMR, personal communication). This could indicate that red king crabs have reached their carrying capacity into the fjord.

Porsangerfjorden, as other fjords in Norway and Svalbard, has experienced significant changes in the last 40 years, with the collapse of the spring spawning herring, cod, and saith fishery, the burst in sea urchin that down grazed macroalgae in the 1970s, the arrival of red king crabs in the 2000s, the recent recovery of macroalgae (Strand, 2019) and the effects of climate change becoming tangible. In this study, we showed that fish assemblages changed with climate change, with an increase of the warm water fishes (significant increase in biomass of Sebastes and increase in the presence frequency of argentine Argentina, blue whiting Micromesistius poutassou, lanternfishes Myctophiformes and round ray Rajella fyllae) at the detriment of cold-water fish (significant decrease in biomass of haddock Melanogrammus aeglefinus and Atlantic hookhead Artediellus atlanticus). The changes undergone by megabenthic assemblages, on the other hand, were caused by the ongoing invasion of red king crabs that lead to a decrease in the abundance of megabenthic taxa and in particular brittle stars Ophiuridae and Iceland scallops Chlamys. The cold-water megabenthic assemblage of the inner part of the fjord should be monitored with caution as it can be closed to the tipping point where red king crab had overgrazed the whole assemblages and durably disrupt their functioning. Overall, this study suggests that invasive species, before climate change effects, would be a major threat to megabenthic assemblages of the subarctic and arctic fjords.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

LM, LJ, MS, and FV designed the study and edited the manuscript. LJ collected the benthic invertebrate’s data. MS provided information and codes on how to extract NorKyst800 model outputs. LM performed the statistical analysis and wrote the first draft. All authors contributed to the article and approved the submitted version.

Funding

This study was conducted in the frame of the project FACE-IT (The Future of Arctic Coastal Ecosystems—Identifying Transitions in Fjord Systems and Adjacent Coastal Areas). FACE-IT has received funding from the European Union’s Horizon 2020 Research and Innovation Program under (grant agreement no. 869154).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank all who made possible the samplings during the coastal surveys and the benthic surveys in the Porsangerfjord. We are grateful to Berengere Husson for her comments on this work. We also thank Jon Albretsen for sharing the river runoff data, Megan O’Sadnick for her insights into sea ice formation into the Porsanger, Ingrid Johnsen Askeland for sharing her data on the water exchanges periodicity between the fjord and open ocean, Stig Dalsøren for his relevant comments onto the NorKyst800 model outputs evaluation, Hans Kristian Strand and Torstein Pedersen for sharing their knowledge of the fjord.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.822979/full#supplementary-material

Footnotes

1. ^ https://www.gebco.net/data_and_products/gridded_bathymetry_data/arctic_ocean/
2. ^ http://geo.ngu.no/download/ShoppingServlet

References