A three-dimensional view on biodiversity changes: spatial, temporal, and functional perspectives on fish communities in the Baltic Sea

Abstract

Fisheries and marine ecosystem-based management requires a holistic understanding of the dynamics of fish communities and their responses to changes in environmental conditions. Environmental conditions can simultaneously shape the spatial distribution and the temporal dynamics of a population, which together can trigger changes in the functional structure of communities. Here, we developed a comprehensive framework based on complementary multivariate statistical methodologies to simultaneously investigate the effects of environmental conditions on the spatial, temporal and functional dynamics of species assemblages. The framework is tested using survey data collected during more than 4000 fisheries hauls over the Baltic Sea between 2001 and 2016. The approach revealed the Baltic fish community to be structured into three sub-assemblages along a strong and temporally stable salinity gradient decreasing from West to the East. Additionally, we highlight a mismatch between species and functional richness associated with a lower functional redundancy in the Baltic Proper compared with other sub-areas, suggesting an ecosystem more susceptible to external pressures. Based on a large dataset of community data analysed in an innovative and comprehensive way, we could disentangle the effects of environmental changes on the structure of biotic communities—key information for the management and conservation of ecosystems.

Baltic Sea, demersal fish community, functional traits, multivariate analysis, spatio-temporal dynamics, species distribution

Introduction

Understanding the impact of environmental conditions on the dynamics and diversity of fish communities is an essential preliminary step for a better prediction of their responses to future changes (Burrows et al., 2011) and for integrative ecosystem-based management (Pikitch et al., 2004; Levin et al., 2009; Möllmann et al., 2014). However, changing environmental conditions can impact biotic communities in multiple ways, and be responsible for changes in structure and function of ecosystems (McGill et al., 2006; Conversi et al., 2015). Environmental conditions are reported to shape the spatial distribution of species (Perry et al., 2005; Poloczanska et al., 2016; Smolinski and Radtke, 2016), influence the temporal dynamics of communities (Rouyer et al., 2008; Möllmann et al., 2009; Hiddink and Coleby, 2012), and select or favour some functional traits (Brind’Amour et al., 2011; Wesuls et al., 2012; Asefa et al., 2017). To the best of our knowledge, no holistic empirical study has investigated simultaneously the effects of environmental changes on (i) spatial distribution, (ii) temporal dynamics, and (iii) functional structure of species assemblages likely due to a lack of appropriate statistical methodologies.

The development of multivariate statistical analyses during the past 20 years has provided ecologists with tools to comprehensively analyse community data and investigate the link between species assemblages, environmental conditions, and functional traits (Dray et al., 2003; Dray and Dufour, 2007; Legendre and Legendre, 2012). Most notably, two frameworks were developed to extend the multivariate methods traditionally limited to the study of the common structure of a pair of data tables (e.g. matrices of species abundance and environmental data, Supplementary Table S1). First, the pair of data tables was extended to study a sequence of paired tables, sequence that could represent different times or spatial locations (Thioulouse et al., 2004; Thioulouse, 2011). These approaches proved to bring new insights into the spatio-temporal structuring of ecological communities (Mazzocchi et al., 2012; Kidé et al., 2015; Chamaille-Jammes et al., 2016). Second, the pair of data tables was extended to a triplet of data tables (Dray and Legendre, 2008; Pavoine et al., 2011; Dray et al., 2014) and allowed the analysis of additional information about traits to discern the traits selected by environmental conditions (Brind’Amour et al., 2011; Wesuls et al., 2012; Asefa et al., 2017). However, these two frameworks are often used separately, limiting studies to investigate the effect of the environment either on the spatio-temporal dynamics of communities or on the selection of traits.

Here, we developed a comprehensive framework based on complementary multivariate statistical methodologies to simultaneously investigate the effects of environmental conditions on the spatial, temporal and functional dynamics of species assemblages, using the Baltic Sea fish community as a case study. The Baltic Sea is a semi-enclosed sea strongly affected by anthropogenic pressures and climate change (Möllmann et al., 2009; Korpinen et al., 2012; Andersen et al., 2015). A strong west-east salinity gradient (Figure 1a) allows the coexistence of approximatively 200 fish species (Ojaveer et al., 2010) ranging from marine to limnic species (Bonsdorff, 2006). The species assemblages are dominated by clupeids, sprat (Sprattus sprattus) and herring (Clupea harengus), that together with cod (Gadus morhua) and flounder (Platichthys flesus) account on average for 90% of the catches. Furthermore, regional climate change models predict an increase in temperature and a decrease in salinity but a high uncertainty remains about the impact of climate change on fish stocks (MacKenzie et al., 2007; Hiddink and Coleby, 2012; Niiranen et al., 2013). The recent increase of anoxic and hypoxic areas in the central Baltic Sea also creates additional pressure on the demersal fish communities (Hinrichsen et al., 2011; Casini et al., 2016; Neumann et al., 2017). Therefore, it is urgent to investigate the role of environmental condition on the spatio-temporal dynamics and structures of fish assemblages in this area.

Figure 1.

Overview of the study. (a) Map of the Baltic Sea with surface salinity of January 2015 represented in scales of blue. (b) Schematic representation of the different dataset and the multivariate methods used for their analysis. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.

Here we provide a coherent and comprehensive analysis of spatial, temporal, and functional dynamics of an entire fish community using modern multivariate statistical approaches. In our case study, we identified sub-assemblages of the Baltic Sea fish community that co-exist and are susceptible to similar environmental conditions. Comprehensive multivariate statistical analyses as suggested and demonstrated here provide crucial information needed for coherent ecosystem-based management of the oceans.

Material and methods

Fish abundance data

Abundance data were collected during the Baltic International Trawl Survey (BITS) (ICES, 2014). Since 2001, the survey has been carried out with a harmonized sampling scheme and a standard gear to sample the demersal fish community in the Baltic Sea. This sampling scheme consists of trawl hauls with a duration of 30 min on average, carried out at a speed of three knots with a demersal otter trawl gear best suited for sampling demersal fish such as gadoids and flatfish. Because of poor survey coverage in shallow areas, we excluded hauls carried out at depths shallower than 20 m. We only included valid hauls carried out during the first quarter of the year (15 February–31 March). In total, our dataset included information from 4086 hauls carried out between 2001 and 2016. With around 250 hauls per year on average, the sampling has a good spatial coverage of the area defined from the Kattegat to the northern Baltic proper (no hauls were recorded at latitude higher than 59°N). The original dataset was downloaded from the ICES Database for Trawl Surveys (DATRAS; http://datras.ices.dk/Home/Default.aspx; data downloaded on the 22 June 2017).

We conducted pre-processing checks to fix mis-identified species and removed non fish species following previous recommendations (Fung et al., 2012; Hiddink and Coleby, 2012). We aggregated species to genus or family level when species were not consistently identified (only for Gobiidae and Callionymus spp.) as suggested by Fung et al. (2012). Among the pelagic species, only sprat and herring were retained in the analyses since they are the two species with the highest catches in the BITS and are consistently caught in almost all the hauls. Other pelagic fish species are not properly sampled by the otter trawl gear and were removed from the dataset. Details of data cleaning are given in Supplementary Table S2. Furthermore, we excluded rare species that occurred in less than 1% of our dataset (i.e. recorded in <40 hauls). The procedure identified 60 rare species (61% of all recorded species), that in total correspond to <0.01% of the total abundance. The final community dataset consisted of abundance values (expressed in number per trawling hour) of 33 species from 4086 hauls.

Traits dataset

Information on the life history of Baltic fish species was retrieved from a previous study (Pecuchet et al., 2016) and the Fishbase database (Froese and Pauly, 2017). Specifically, trophic level and maximum length (in cm) were extracted from Fishbase (Froese and Pauly, 2017). Mean fecundity (in number of eggs spawned per adult female in one spawning season), the shape of the caudal fin (in five categories: rounded, truncated, emarginated, forked and continuous) and the body shape (in 4 categories: gadoid-like, flat, elongated, and eel-like) were obtained from Pecuchet et al. (2016). In total, two qualitative and three quantitative traits were used to characterize the 33 fish species. These five different traits are complementary and describe the diet, habitat and reproduction for all species included in the study (Pecuchet et al., 2016). For graphical visualization, species were identified in accordance with their FAO 3-Alpha Species Codes provided in Table 1 (http://www.fao.org/fishery/collection/asfis/en, version February 2017).

Table 1.

Fish species found in the Baltic Sea, ordered per spatial cluster.

	Species	FAO name	Trophic level	Maximum length (cm)	Fecundity (no. eggs)	Caudal shape	Body shape
Kattegat	Arnoglossus laterna	MSF	3.6	25	50 000	Rounded	Flat
	Callionymus spp.	YVX	3.3	32	5000	Rounded	Gadoid-like
	Eutrigla gurnardus	GUG	3.6	60	245 000	Emarginated	Gadoid-like
	Glyptocephalus cynoglossus	WIT	3.1	60	100 000	Rounded	Flat
	Hippoglossoides platessoides	PLA	3.5	83	380 000	Rounded	Flat
	Lepidorhombus whiffiagonis	MEG	4.2	60	500 000	Rounded	Flat
	Lumpenus lampretaeformis	LMJ	3.6	50	700	Rounded	Eel-like
	Melanogrammus aeglefinus	HAD	4	112	9 000 000	Emarginated	Gadoid-like
	Merluccius merluccius	HKE	4.4	140	1 000 000	Truncated	Elongated
	Microstomus kitt	LEM	3.2	65	200 000	Rounded	Flat
	Pholis gunnellus	FGN	3.5	25	100	Rounded	Eel-like
	Scophthalmus rhombus	BLL	3.8	75	5 000 000	Rounded	Flat
	Solea solea	SOL	3.1	70	300 000	Rounded	Flat
	Trachinus draco	WEG	4.2	53	57 600	Truncated	Gadoid-like
	Trisopterus esmarkii	NOP	3.2	35	220 000	Emarginated	Gadoid-like
	Trisopterus minutus	POD	3.7	40	10 000	Truncated	Gadoid-like
Western	Agonus cataphractus	AFT	3.4	21	3000	Rounded	Elongated
	Gobiidae spp.	FGX	3.2	10	3000	Rounded	Elongated
	Limanda limanda	DAB	3.3	40	150 000	Rounded	Flat
	Merlangius merlangus	WHG	4.2	70	400 000	Truncated	Gadoid-like
	Pleuronectes platessa	PLE	3.3	100	552 000	Rounded	Flat
	Pollachius virens	POK	4.2	130	2 900 000	Emarginated	Gadoid-like
Eastern	C. harengus	HER	3.2	45	60 000	Forked	Gadoid-like
	Cyclopterus lumpus	LUM	3.8	61	100 000	Truncated	Gadoid-like
	Enchelyopus cimbrius	ENC	3.5	41	500 000	Rounded	Elongated
	Gadus morhua	COD	4.3	200	1 000 000	Truncated	Gadoid-like
	Gasterosteus aculeatus	GTA	3.4	11	350	Truncated	Gadoid-like
	Myoxocephalus quadricornis	TGQ	3.7	60	18 000	Truncated	Gadoid-like
	Myoxocephalus scorpius	MXV	3.6	60	10 000	Truncated	Gadoid-like
	Platichthys flesus	FLE	3.2	60	1 000 000	Rounded	Flat
	Scophthalmus maximus	TUR	4	100	5 000 000	Rounded	Flat
	Sprattus sprattus	SPR	3	16	10 000	Forked	Gadoid-like
	Zoarces viviparus	ELP	3.5	52	100	Continuous	Elongated

Traits are derived from Pecuchet et al. (2016, 2017) and Froese and Pauly (2017).

Environmental dataset

Based on the location and time of each of the 4086 hauls, we extracted nine environmental variables. We selected environmental variables based on their potential or known effect on the demersal fish community, specifically depth, local hydrographic conditions, primary productivity and large-scale climatic conditions. Additionally, we assured that the selected variables were not strongly cross-correlated (Pearson coefficient < 0.7).

Trawl depth was retrieved directly from the information provided in the BITS dataset. Five variables concerning local hydrographic conditions were derived from the Baltic Sea Ice-Ocean Model (BSIOM) (Lehmann and Hinrichsen, 2000; Lehmann et al., 2002, 2014), a hydrodynamic model with an oxygen consumption calculation sub-module. BSIOM provided values of temperature, salinity, and oxygen with a horizontal resolution of 2.5 km and 60 vertical levels. The temporal evolution of three-dimensional temperature, salinity, and oxygen fields are in good agreement with hydrographic measurements of the ICES database (Lehmann et al., 2014). The five variables used in this study were annual bottom temperature, oxygen and salinity, surface temperature at the time of the survey and seasonality of bottom temperature. Annual bottom hydrographic conditions were selected as characteristic of the habitat (temperature: sbt_an, salinity: sbs_an, and oxygen: oxb_an). Average surface temperature during the first quarter (sst_q1) was selected as a snapshot of the hydrographic conditions at the time of the survey. The seasonality of bottom temperature (sbt_ra) was estimated as the range between average monthly temperatures and is an indicator of seasonal variation of the benthic habitat. Two variables concerning primary productivity were estimated from the Chlorophyll a concentration (mg.m⁻³) of the GlobColour project (Maritorena et al., 2010), merging Ocean Colour products from different sensors. We used a monthly averaged dataset with a spatial resolution of 1 km, downloaded from http://hermes.acri.fr/on 13 June 2017. The two variables selected as indicators of the primary production were Chlorophyll a concentrations averaged over the first quarter (chl_q1) and over the previous year (chl_an). Large-scale climate conditions were represented in our analysis by the North Atlantic Oscillation (NAO) index (Hurrell, 1995), which indicates high frequency (7–25 years) atmospheric variations and is known to affect Baltic biotic communities (Hänninen et al., 2000; Möllmann et al., 2009). The NAO time series was downloaded from the climate indices platform of the Earth System Research Laboratory: http://www.esrl.noaa.gov/psd/data/climateindices/list/ on 2 June 2017.

Multi-tables multivariate analyses

First, each dataset (i.e. species abundance, species traits and environmental data) was analysed according to appropriate ordination methods corresponding to the nature of the variables (Figure 1b, Supplementary Table S1). Correspondence analysis (CA) was computed on the fish abundance dataset (a matrix of species by sample). CA is well suited for abundance data along large environmental gradient because species communities often show a unimodal distribution along a gradient and, using the chi-square distance, CA can highlight differences of species composition profiles (Legendre and Gallagher, 2001; Greenacre, 2017). Abundance was previously log-transformed (x + 1) to reduce the influence of the dominant species in the analysis of community structure. Principal component analysis (PCA) was performed on the environmental dataset with nine quantitative variables (a matrix of environment variables by sample) using the row weights (corresponding to the samples) from the previous CA on the fish abundance dataset in order to permit the comparison between species distribution and environmental conditions. The trait dataset (a matrix of species by trait) contained a mix of quantitative and qualitative variables and was analysed using the Hill and Smith (1976) method. This method combined a PCA on quantitative variables and a multiple CA on qualitative variables. Species were weighted according to column weights in the previous CA on the fish abundance dataset, in order to permit the comparison between species distribution and species traits.

Second, we used within and between-group analysis to assess and separate spatial and temporal variabilities, with the year of sampling as grouping variable (Dolédec and Chessel, 1987; Franquet et al., 1995). Between-group analysis is analogous to an ordination of the table of group means. In other words, it seeks to reveal the main temporal pattern by looking for the highest differences among years. Within-group analysis is the reverse of between-group analysis, i.e. it is the ordination of the residuals among initial data and group means. It removes the effect of the grouping variable and analyses the remaining variability, so in our case the spatial variability.

Third, co-inertia analysis (COA) was used to link fish community composition with environmental conditions by coupling these two data tables (Dray et al., 2003, 2014). COA is an unconstrained symmetric analysis that searches for axes that maximize the covariance between the samples of both data tables. We applied COA on the results of the between and within-group analysis, resulting in the so called Between-group co-inertia analysis (BGCOA) and the Within-group co-inertia analysis (WGCOA) (Thioulouse, 2011). BGCOA reveals the temporal co-dynamics between fish species and the environment. WGCOA shows the spatial structure of the fish community that could be explained by environmental conditions. Beforehand, the association between the two tables was tested using a Monte-Carlo permutation test and the RV coefficient (Heo and Gabriel, 1998). The RV coefficient is a generalization of the Pearson’s correlation coefficient for matrices (instead of vectors). A permutation test with 1000 random permutations was performed to evaluate if the association between the two data tables was significantly stronger than expected by chance. In the cases that the p-value was higher than 0.05, the results of the COA were not analysed and are not presented in the “Results’ section.

Finally, fourth-corner and RLQ methods were used to assess the link between the species trait composition and environmental variation (Dolédec et al., 1996; Dray and Legendre, 2008; Dray et al., 2014). The fourth-corner method is a permutation test to evaluate the pairwise association between traits and environmental variables, measured by a Pearson’s correlation coefficient. We used a combination of permutations of samples and of species to correct for inflated type I error (Dray and Legendre, 2008). RLQ is a multivariate method that assesses the trait-environment relationships (Dolédec et al., 1996). Partial RLQ takes into account the partition of environmental variation in within and between-groups analysis (Wesuls et al., 2012). In the case that none of the relationships between species traits and the environment was significant in the fourth-corner test, the results of the RLQ analysis are not presented in the Results section.

Definition and characterization of spatial sub-assemblages

We defined sub-assemblages of Baltic Sea fish species that share similar spatial distributions and hence are favoured by similar environmental conditions. We computed Euclidean distances between fish species from the projection of species on the PCs of the WGCOA, and subsequently conducted a hierarchical cluster analysis based on Ward’s criterion (Ward, 1963). Based on a graphical interpretation of the dendrogram, we selected the number of clusters. The robustness of the selected number of clusters and of the clustering solution was tested by comparison with the alternative k-means cluster analysis. The clustering provided a simplification of the fish community into fewer sub-assemblages sharing similar spatial distributions. We characterize these sub-assemblages by looking at their temporal dynamics and functional richness. Functional richness was calculated as the area of the convex hull on the functional space, i.e. the volume of the functional space occupied by the community (Villéger et al., 2008). The functional space is defined from species projections on the principal components of the previous Hill and Smith analysis.

Software and sources

All statistical analyses were conducted in the programming environment R 3.3 (R Core team, 2017). The ade4 package (Dray and Dufour, 2007) was used to compute the multivariate analyses. The functional richness was calculated with the FD R package (Laliberté and Legendre, 2010). Maps were created with the mapdata package (Becker et al., 2016). The cleaned datasets and the R-script are available in Supplementary Material S1.

Results

Spatial distribution of the fish community linked with environment

The structure of the fish community in the Baltic Sea was strongly linked with salinity conditions and depth (Figure 2). The first two principal components of the within-group co-inertia analysis (WGCOA) explained 95% of the covariance between fish abundance and environmental variables. The first principal component (PC1, 87% of the covariance) separated fish species favouring highly saline waters in the Kattegat, against fish species inhabiting less saline waters in the Baltic Proper (Figure 2a and b). The salinity gradient was also associated with higher bottom temperatures (sbt_an), higher primary production (chl_q1) and shallower waters (Figure 2f). Most of the fish species had a negative score on PC1 (Figure 2e, left side on the x-axis), i.e. were located in highly saline waters. The second PC, explaining only 8% of the covariance, represented mainly the differences between shallow and deep waters (Figure 2c and d). Deep basins were also associated with lower seasonal variation in bottom temperature and lower oxygen content (Figure 2f). Some species strongly preferred shallow waters (Scophtalmus maximus, TUR), while others were caught mainly in deep basins (Enchelyopus cimbrius, ENC) (Figure 2e).

Figure 2.

The within-group co-inertia analysis summarized the spatial links between community composition and environmental conditions in two PCs. The red to blue colour gradient represents the hauls’ score on the PC1 based on fish composition (a) or environmental variables (b). The purple to green gradient represents the hauls’ score on PC2 based on fish composition (c) or environmental variables (d). The link between fish composition and environmental conditions can be visualized by the scores on the two first PCs of fish species (e) and environmental variables (f). Species are represented following the three-letters code shown in Table 1. The names of the environmental variables showing low scores on the PCs were abbreviated: sbt_an is the annual bottom temperature, sst_q1 is the surface temperature during the first quarter, chl_q1 and chl_an are chlorophyll a concentrations averaged over the first quarter and over the previous year, respectively.

Three sub-assemblages of the Baltic Sea fish community were identified with hierarchical clustering analysis based on their PC scores derived by the WGCOA (Figure 3a, as also confirmed by k-means, Supplementary Figure S1). The sub-assemblages grouped species according to their spatial distribution and to whether they were sharing similar environmental conditions. The differences between sub-assemblages were mainly defined along PC1 (Figure 3b, x-axis), in other words were strongly linked to the west-east salinity gradient. A sub-assemblage of 16 Kattegat fish species (Table 1) were favoured by high saline waters, therefore inhabiting only the Kattegat (latitude higher than 56°N and longitude lower than 13°E). The Kattegat sub-assemblage included among the most abundant species, American plaice (Hippoglossoides platessoides, PLA), Norway pout (Trisopterus esmarkii, NOP), dragonets (Callionymus spp., YVX), and greater weever (Trachinus draco, WEG). The cluster analysis identified another sub-assemblage of 6 western Baltic fish species (Table 1), adapted to middle salinity conditions, with a distribution ranging from the Kattegat to the Arkona basin (longitude lower than 15°E). The western Baltic fish species included dab (Limanda limanda, DAB), whiting (Merlangius merlangus, WHG) and European plaice (Pleuronectes platessa, PLE). The third sub-assemblage comprised 11 eastern Baltic fish species (Table 1) that were favoured by low salinity conditions and could potentially inhabit the entire study area, from the Kattegat to the northern Baltic proper. The eastern Baltic group included sprat (S. sprattus, SPR), herring (C. harengus, HER), cod (G. morhua, COD) and flounder (P. flesus, FLE).

Figure 3.

Three sub-assemblages were identified from a cluster analysis of Baltic Sea fish species according to their spatial distribution. (a) Dendrogram of the cluster analysis, suggesting three distinct groups. (b) Projection of the groups on the two first PC of the within-group co-inertia analysis. The key environmental drivers are shown on the PC. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.

Temporal and functional patterns of the fish community

According to the between-group analysis, the temporal dynamics of the Baltic Sea fish community accounted for only 2% of the variance of fish abundances, suggesting a relatively stable structure of fish assemblages during the period 2001–2016. The relationship with environmental dynamics, tested with a permutation test using the RV coefficient, was not significant (p-value = 0.1) (Supplementary Figure S2). Therefore, the results of the between-group co-inertia analysis linking fish dynamics and environment are not presented here, but rather the results of the between-group analysis of fish community dynamics (even if representing only 2% of total spatio-temporal variance). The main mode of variability in the fish community dynamics was associated with a general increase in species abundances between 2001 and 2016 (PC1, explaining 39% of the temporal variance) (Figure 4a). The species with the highest relative increase were Arnoglossus laterna (MSF) and Myoxocephalus quadricornis (TGQ) (Figure 4b). Some species also experienced a decrease, especially Lepidorhombus whiffiagonis (MEG) that was last recorded in the Baltic Sea in 2011. The second PC, explaining 19% of the temporal variance, highlighted the difference in abundances between the years 2001 and 2013 and the year 2016, which was mainly characterizing the dynamics of saithe (Pollachius virens, POK) (Figures 4a and b).

Figure 4.

Temporal and functional characteristics of fish community. (a) The temporal dynamics revealed by betweengroup analysis and summarized in 2 PCs, explaining respectively 39 and 19% of the total variance. (b) Projection of fish species on these temporal PC is displayed with colour associated to their spatial subassemblage. (c) Functional space, revealed by Hill and Smith analysis, with 2 PCs explaining, respectively 32 and 23% of the total variance. (d) Projection of fish species on these functional PC is displayed with colour associated to their spatial sub-assemblage. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.

The ordination of the five functional traits with the Hill and Smith analysis highlighted two main functional characteristics of fish species along the first two PCs, together explaining 55% of the traits variance (Figure 4c). The first PC (32% of the total variance) summarized the small-large continuum. Large, high trophic level, and high fecundity species such as cod (G. morhua, COD) were separated from species that are small, with low trophic level and low fecundity such as rock gunnel (Pholis gunnellus, FGN) or sprat (S. sprattus, SPR) (Figure 4d). The second PC (23% of the total variance) revealed the difference between flat fish species with rounded caudal fin and the gadoid-like shaped species. The link between fish traits and environmental variables was tested with a fourth corner permutation test and no significant relation was found between individual traits and environmental variables.

Sub-assemblages characteristics and local biodiversity indices

The spatial distributions of the three sub-assemblages were nested, i.e. the western Baltic sub-assemblage also inhabits the Kattegat, and the eastern Baltic sub-community was present all over the surveyed area (Figure 5a). The temporal dynamics were quite diverse but, on average, the abundance of fish species had increased in the observation period (grey shaded area in Figure 5b). However, we observed differences between the sub-assemblages, the Kattegat displaying the lowest relative increase in abundances (apart of 2016) and the eastern Baltic sub-assemblage the highest. Interestingly, the functional richness of the Kattegat and eastern Baltic sub-assemblages was high (Figure 5c).

Figure 5.

Spatial, temporal and functional characterization of the sub-assemblages. (a) Average spatial distribution, with the intensity of colours proportional to the spatial distribution. (b) Temporal dynamics from 2001 to 2016. The bold line represents the median relative abundance, the shaded area the inter-decile range. The whole fish assemblage is represented in grey, the sub-assemblages in their respective colours. (c) Functional richness of the sub-assemblages, compared with the whole community; PC1 in x-axis represents the difference between large (left) and small (right) fish species, PC2 in y-axis represent the difference .between flat (up) and gadoid-like shaped (bottom) fish species. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.

The spatial overlap of the three sub-assemblages was confirmed by the community composition per haul (Figure 6a and b). Hauls in the Kattegat (defined by latitude higher than 56°N and longitude lower than 13°E) were composed, on average, of 34% of Kattegat fish, 31% of western Baltic fish and 35% of eastern Baltic fish. Hauls carried out in the Baltic Proper (longitude >15°E) were nearly exclusively composed of species from the eastern Baltic sub-assemblage (95 and 5% from the western Baltic sub-assemblage). The spatial distribution of species richness per haul confirmed the increase of species richness along the salinity gradient (Figure 6c). As expected from the nested sub-assemblages and the high functional richness of the eastern Baltic sub-assemblage, functional richness had a relatively lower variation along the salinity gradient (Figure 6d). The recent increase in abundance of some fish species of the eastern Baltic assemblage could explain the recent increase of species richness among the less diverse hauls (the ninth decile showed an increase since 2009, while the median is more or less stable) (Figure 6e).

Figure 6.

Fish biodiversity in the Baltic Sea from the information of more than 4000 samplings. (a, b) Hauls composition of the three sub-assemblages are represented in RGB colour scale with red for the Kattegat, green for the Western and blue for the Eastern sub assemblages. (c, d) Spatial distribution of species and functional richness per haul. (e) Temporal evolution of species richness, the line represents the median, the shaded area represents the interdecile range.

Discussion

Environmental conditions drive fish community composition

Based on a large dataset of more than 4000 samples and using complementary multivariate analyses and statistical tests, we investigate the links between fish communities and environmental conditions in the Baltic Sea. Salinity, decreasing from marine waters in the Kattegat to brackish waters in the Baltic Proper, is the main driver of fish community composition. Along this gradient, our statistical approach is able to identify three sub-assemblages within the overall Baltic fish community. The three sub-assemblages are nested, with most fish species inhabiting the Kattegat. This finding agrees with the predicted reduction of species richness from marine to brackish water (Pecuchet et al., 2016; Smolinski and Radtke, 2016), also found in the Baltic Sea zoobenthic community (Bonsdorff, 2006). Additionally, our study disentangles the species composition and identifies 16 species strongly limited in their spatial distribution by salinity, preferring high-salinity conditions (the Kattegat sub-assemblage). Interestingly, the depth gradient is often reported as the most important environmental driver shaping the fish community in other large marine ecosystem (Kidé et al., 2015; Dencker et al., 2017; Pecuchet et al., 2017). Here we find a weaker linkage of fish assemblages with depth, confirming the very unique conditions in the semi-enclosed brackish waters of the Baltic Sea. The salinity gradient is, by far, the main driver of fish assemblages, suggesting that the Baltic Sea could be more similar to a large estuary than open ocean. Although this information is not novel, our study compares both drivers quantitatively and we find that salinity explain 87% of the covariance between fish and environmental conditions, while the depth gradient accounts only for 8%. If similar methodological framework would be applied to other ecosystems, we could compare the importance of different drivers across marine ecosystems. Our approach needs a large amount of collected data, which are already available for intensively monitored seas in Europe (Granger et al., 2015; Dencker et al., 2017; Frainer et al., 2017) or in North America (Batt et al., 2017).

Mismatch between taxonomic and functional diversity

Even though species distributions are highly linked with environmental gradients, we do not find any significant relationship between the functional characteristics of fish species and environmental conditions. If species would be selected randomly, the functional richness would tend to increase with the number of species (Mouillot et al., 2007). On the contrary, the spatial overlap of sub-assemblages and the high functional richness of Kattegat and eastern Baltic sub-assemblages suggest that the number of species is reduced along the west-east gradient but without a decrease of functional richness (Figures 5c and 6c). This is especially surprising for the eastern Baltic sub-assemblage, which includes the few species that can tolerate the low salinity conditions. These remaining species are able to occupy all the “niches” defined in the functional space, suggesting that the environmental conditions may limit similarities between the remaining species thus favouring the realisation of all the niches needed for the functioning of communities. This result agrees with Pecuchet et al. (2016) that proved a distinction between environmental filtering acting in the western Baltic Sea and neutral or limiting similarity acting in the Baltic Proper. However, Pecuchet et al. (2016) also found a link between functional richness and salinity when the diversity indicator was aggregated spatially into a regular grid. This link is not confirmed by our analysis made at the species level and considering each individual haul. The difference can be explained by some outlier hauls in the Baltic Proper with low catches, resulting in an abnormal low functional richness that can have a high influence on spatially averaged values. Moreover, the limited number of traits, although usual in functional studies of fish assemblages (Dencker et al., 2017; Pecuchet et al., 2017), covers only the life history strategies (survival, growth, reproduction) and do not take into account tolerance range of species (e.g. temperature or salinity preferences). Adding environmental tolerance traits would clearly increase the link between environment and traits, but our goal was to focus only on the life history strategies.

Limits of data-driven approach

As with any statistical approach, the ability of the methods applied here is limited by the quality and amount of data available. For example, the dataset used covers the 16-year period from 2001 to 2016, i.e. after the regime shift occurring in the Baltic Sea during the late 1980s (Möllmann et al., 2009; Casini, 2013). Including the period prior to the shift would likely increase the importance of the temporal dynamics and the capacity to detect a significant link with environmental variability. Even though the sampling started in 1991, we did not include data prior to 2001 because the sampling was performed using different gears and the sampling scheme of the surveys was different, potentially affecting the robustness of our analysis. The short length of the time series is the main limitation of our study, and it stresses the importance of rigorous and continuous data collection, a very valuable source of information in order to understand, preserve and manage marine ecosystems in a better way. In our analysis, the absence of a link between the temporal dynamics of fish species and the environmental conditions in the period 2001–2016 is informative. This finding is contrary to Hiddink and Coleby (2012) that linked the dynamics of species richness with temperature in Kattegat and salinity in Baltic Proper. The difference can be explained by a different time period (1990–2008) used by Hiddink and Coleby (2012) and the fact that our approach does not aggregate species into a diversity indicator and assume homogeneous dynamics over the whole study area. While looking at the dynamics of the sub-assemblages, we find that Kattegat and eastern Baltic sub-assemblages have different dynamics, suggesting the use of methods that could study the interaction between spatial distribution and temporal dynamics, such as tensor decomposition (Frelat et al., 2017).

Moreover, we could not include fishing pressure or other direct anthropogenic pressures in our study because they are difficult to estimate and not available at the spatial resolution of our analysis. Yet, it is clear that human pressures have a strong impact on the fish community (Korpinen et al., 2012; Andersen et al., 2015) and their trait composition (D’agata et al., 2014; Henriques et al., 2014; Koutsidi et al., 2016). Another limit of our study remains in the fact that we study abundance at species level, which may hide information about different sub-populations of the same species. For example, cod is known to be divided in two populations: the eastern and western Baltic cod stocks (Aro, 1989; Bagge et al., 1994), but here included as only one species, which may blur the temporal dynamics of these two stocks.

Management implications and concluding remarks

The mismatch between taxonomic and functional diversity, associated with the spatial overlap of sub-assemblages, suggests that the functional redundancy decreases from west to east in our study area. The low functional redundancy in the Baltic Proper implies that its ecosystem is susceptible to changes in external pressures such as hydrography, nutrient inputs, and fisheries overexploitation that can provoke drastic reductions in fish abundances (Rice et al., 2013). Therefore, fisheries management in the Baltic Proper should be precautious by taking in consideration the specific local characteristics of the fish community. Our study demonstrates that based on a large dataset of community data, analysed in an innovative and comprehensive way, we can provide a complete view of the effects of environment on the structuring of biotic communities in space, time and functions. Similar methodological framework can be used in other large marine ecosystems to gain better understanding of the effect of environmental variations on biodiversity, key information for the management and conservation of ecosystems.

Acknowledgements

The authors are grateful to all contributors of the DATRAS database who collected and took part in the Baltic International Trawl Survey, and to ICES who merged, cleaned and made the data open access. We thank the two anonymous referees and Marta Coll for their useful comments that helped improve this manuscript. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 675997 (ITN MARmaED). The results of this publication reflect only the author's view and the Commission is not responsible for any use that may be made of the information it contains. The study is a contribution to the University of Hamburg Center of Excellence on “Integrated Climate System Analysis and Prediction” (CliSAP) funded by the German Science Foundation. This work was also financed by the BONUS INSPIRE project supported by the Joint Baltic Sea Research and Development Programme BONUS (Art 185), funded jointly by the EU and the Swedish Research Council Formas.