Habitat Heterogeneity Determines Climate Impact on Zooplankton Community Structure and Dynamics

Understanding and predicting species distribution in space and time and consequently community structure and dynamics is an important issue in ecology, and particularly in climate change research. A crucial factor determining the composition and dynamics of animal populations is habitat heterogeneity, i.e., the number of structural elements in a given locality. In the marine pelagic environment habitat heterogeneity is represented by the distribution of physical oceanographic parameters such as temperature, salinity and oxygen that are closely linked to atmospheric conditions. Little attention has been given, however, to the role of habitat heterogeneity in modulating the response of animal communities to external climate forcing. Here we investigate the long-term dynamics of Acartia spp., Temora longicornis, and Pseudocalanus acuspes, three dominant zooplankton species inhabiting different pelagic habitats in the Central Baltic Sea (CBS). We use the three copepods as indicator species for changes in the CBS zooplankton community and apply non-linear statistical modeling techniques to compare spatial population trends and to identify their drivers. We demonstrate that effects of climate variability and change depend strongly on species-specific habitat utilization, being more direct and pronounced at the upper water layer. We propose that the differential functional response to climate-related drivers in relation to strong habitat segregation is due to alterations of the species’ environmental niches. We stress the importance of understanding how anticipated climate change will affect ecological niches and habitats in order to project spatio-temporal changes in species abundance and distribution.


Introduction
Understanding and predicting species distribution in space and time and consequently community structure and dynamics is an important issue in ecology, and particularly in climate change research. Animal populations in terrestrial and aquatic environments shift their geographical range or show strong fluctuations in abundance in response to climate-induced changes in, e.g., temperature, precipitation or ocean circulation [1][2][3][4][5]. Changes on the individual species and population level can eventually cause community re-organizations on large spatial scales [6][7][8][9] with concomitant effects on ecosystem functioning [10].
A crucial factor determining the composition and dynamics of animal populations is habitat heterogeneity, i.e., the number of structural elements in a given locality [11,12]. For instance, species diversity is known to increase with habitat heterogeneity [13][14][15][16], since structurally diverse habitats promote the coexistence of a higher number of species by providing more niches and ways of exploiting environmental resources [17,18]. In aquatic systems, mostly submerged plants and substrate material determine the physical structure of the environment, and therefore, have a considerable influence on the distributions and interactions of demersal fish species [14,19], seagrass, stream or sub-and intertidal macroinvertebrates [20][21][22][23]. In the marine pelagic environment however, habitat heterogeneity is represented by the distribution of physical oceanographic parameters such as temperature, salinity and oxygen. Although these hydro-climatic variables are closely tied to atmospheric conditions, the question of how habitat heterogeneity modulates the response of animal communities to external climate forcing has received little attention so far.
The pelagic of the Central Baltic Sea (CBS) is characterized by spatial hydrographic gradients that create a number of different habitats for local animal populations. Horizontally, both salinity and temperature generally decline from West to North-East (Fig. 1C). Even more pronounced, vertical habitat heterogeneity is induced by stratification of the water column during most of the productive season into (i) a warm but low saline surface layer, (ii) a cold and low saline intermediate layer, and (iii) a cold and highly saline deep layer (Fig. 1B). Here we investigate the importance of the spatial habitat structure for the Baltic Sea zooplankton community. Zooplankton is a suitable indicator of the effect of climate change on aquatic ecosystems and long-term changes in marine zooplankton are well studied, e.g. [24][25][26].
The CBS zooplankton community is composed of six major taxonomic groups (copepods, cladocerans, rotifers, appendicular-ian, ctenophores, polychaetes, and bivalves), dominated by very few species (often 1 to 3) [27,28]. Physiological adaptations to physical oceanographic conditions result in different utilization of the different water masses or habitats available [27,28]. In our study we investigate the long-term dynamics of Acartia spp., Temora longicornis, and Pseudocalanus acuspes, three dominant zooplankton species that make up to 60% of the spring biomass in the CBS and distribute differently according to these habitats. Acartia spp. mainly inhabits the surface water, T. longicornis the mid-water and P. acuspes the deeper, saline layer (Fig. 1B) [29,30]. The vulnerability of our target species to predation is as well strongly habitat-specific. The dominant planktivore predators herring (Clupea harengus) and sprat (Sprattus sprattus) feed mainly in deeper waters [31], hence are displaying a different overlap with their zooplanktonic prey species (Fig. 1B).
Here we use the three copepods as indicator species for changes in the CBS zooplankton community, and specifically test (i) whether our target zooplankton species display different population trajectories in space and time, and (ii) if hydro-climatic factors and predation pressure differentially affect the long-term development of the individual species populations. We apply statistical modeling techniques to compare species-, area-, and seasonspecific long-term population trends and to identify their main drivers.

Zooplankton data
For our study we have, for the first time, combined data sets from three long-term zooplankton monitoring programs in the Central Baltic Sea (CBS; Fig. 1A). With this new and unique data set we were able to investigate the CBS zooplankton dynamics within all three major basins of the CBS (the Bornholm Basin (BB), Gdansk Deep (GD), and Gotland Basin (GB)) for the recent five decades . The data sets were provided by the Latvian Institute of Food Safety, Animal Health and Environment (BIOR), the Finnish Environment Institute (SYKE), and the Leibniz-Institute for Baltic Sea Research (IOW) in Germany. The programs differ in gear type and spatio-temporal coverage, with the Latvian program being the most comprehensive and described in detail by Möllmann et al. [25]. For an overview of sampling programs see Table S1.
In order to construct basin-specific time series on spring (sampling mainly in May) and summer (sampling mainly in August) biomass (mg*m 23 ) for Acartia spp., Temora longicornis and Pseudocalanus acuspses, we merged the three data sets according to the following procedure: N We started by calculating the annual mean basin-and seasonspecific biomass of each species for each data set.
N After log10(x+0.001) transformation, we combined the three data sets by calculating annual means for each species, basin and season weighted by the sampling number of each individual sampling program.
N Subsequently, we accounted for differences between sampling programs due to gear-specific capture efficiencies and different spatial and temporal resolutions following the approach suggested by Mackas and Beaugrand [24] and computed species-specific log10(x+0.001) transformed anomalies: where BM sb (y) is the mean biomass of a species in season s and basin b in a particular year y, BM sb represents the season-and basin-specific overall time series mean for each species.
Using log-transformed anomalies accounts for any potential sampling bias that is inherent in each annual mean value as well as in the overall time series mean [24]. It should, however, be noted that this approach does not allow for a quantitative comparison between basins, seasons, or species, but rather provides a basis for comparisons of long-term trends.

Environmental covariates
We included several hydro-climatic variables and predation pressure as covariates in our statistical analyses. The main hydrographic variables affecting Acartia spp., T. longicornis and P. acuspes are water temperature and salinity [25,[32][33][34]. We computed therefore basin-specific time series of seasonal (spring -April/May; summer -July/August) temperature and salinity conditions (Fig. S1). For both time periods, we calculated mean values of individual depth ranges inhabited by each copepod species based on [29,30,33,35], including 0 to 20 m and 20 to 40 m for Acartia spp. and T. longicornis respectively (see Table S2). P. acuspes displays a strong ontogenetic shift in depth distribution (e.g., nauplii and copepodites C1-3 dwelling in the upper and midwaters; copepodites C4-5 and adults dwelling close to or in the deep permanent halocline) and characteristic stage-distributions in each season [33,36]. Hence, we considered three different depth ranges for spring, i.e., 20-60 m for early copepodites and 60-90 m (BB/GD) or 70-100 m (GB) for older stages, depending on the local depth of the permanent halocline. Eventually, we assumed 30-50 m to be the distribution range of the P. acuspes population in summer. Hydrographic data were retrieved from the International Council for the Exploration of the Sea (ICES) website (http:// www.ices.dk/ocean/aspx/HydChem/HydChem.aspx).
As an indicator for the overall climate condition we used the Baltic Sea Index (BSI) provided by the GEOMAR Helmholtz Centre for Ocean Research Kiel, Germany (Table S2, Fig. S2A). This local climate mode represents the difference of normalized sea level pressure anomalies between Oslo in Norway and Szcecin in Poland. Although closely correlated to the North Atlantic Oscillation (NAO), the BSI generally has been better related to local oceanographic processes in the Baltic than remote, largescale atmospheric patterns such as the NAO [37]. We used the winter BSI (December-March) since winter patterns of local and regional wind fields have a strong influence on the hydrographic situation in early spring [38], affecting phytoplankton dynamics and the onset of the spring bloom [39].
Sprat (Sprattus sprattus) is the dominant planktivore in the Baltic Sea basins and, at times, appears capable of exerting top-down control of P. acuspes population [40]. This makes sprat stock size generally a suitable indicator of predation pressure [41]. However, official stock size estimates for sprat are not available prior to 1974 [42]. Hence, a predation index (PI) was developed using stock size estimates for cod (Gadus morhua). Eastern Baltic cod is the major piscivore in the Baltic Sea and its stock size is inversely related to the stock size of sprat [41,43] (Fig. S2B). We used spawning stock biomass (SSB) estimates for cod from official stock estimates that date back to 1966 [42], and extended the time series back to 1960 using estimates by Eero et al. [44]. The combined time series was inverted to mimic the predation pressure by planktivorous sprat and standardized to values between 0 and 1 and. The max. value of cod SSB is represented therefore by a PI value of 0 while the min. cod SSB value is indicated by a PI value of 1.
Food-limitation can be an important influence on population dynamics of copepod species, but studies on the feeding ecology of Baltic copepods in the field are scarce. Some earlier studies indicate a stronger importance of hydro-climatic variables as well as predation for the dynamics of P. acuspes [41,43]. Other recent studies show the importance of food availability for the dynamics of Acartia longiremis. and T. longicornis in different seasons [35,45]. Unfortunately, no reliable estimates of the different food sources were available for our full investigation period. Changes in phytoplankton abundances and timing are often affected by temperature, e.g., through changes in light-saturated photosynthesis rates [46], as well as by irradiance (cloud/solar radiation) and wind-induced vertical mixing intensity [47,48]. Mixing intensity in particular is a key variable for growth performance as it determines largely the light and nutrient conditions [49]. Under westerly wind conditions (represented by positive BSI anomalies), cloudiness is generally higher and wind mixing reduced [37,39]. We therefore assumed that dynamics of food availability are indirectly reflected by the temperature and BSI time series.

Statistical analysis
We used Generalized Additive Modeling (GAM) techniques [50,51] in this analysis which have the advantage of not requiring an a priori specification of the relationship between the response variable (Y) and the explanatory variable (X). Each Y i is here linked with X i by a smoothing function instead of a coefficient b as in traditional regression techniques and hence relationships do not have to be linear.
In a first analysis we tested for spatial differences in temporal trends in species biomass between the three basins. We modeled the log-transformed spring or summer biomass anomalies of each species either as a single function of time (year), i.e., for the entire CBS area, or as basin-specific functions of time where BM s by is the season-specific (index s) log-transformed biomass anomaly in basin b and year y. f and f b are the single and respectively basin-specific thin plate regression spline functions of the variable time. a is the intercept and e by the random noise terms assumed to be normally distributed with zero mean and finite variance. To focus only on the major trends and avoid over-parameterization, we limited the effective degrees of freedom (edf) to a maximum of 3. For each of the species-and seasonspecific comparisons between models including a spatially uniform or basin-specific smoothing functions, we applied an F-ratio test [51]. Significant differences (p-values,0.05) indicated here that the more complex GAMs with basin-specific smoothers performed better and, hence, long-term trends are better described for each basin separately.
In a second analysis we used the non-linear regression technique to identify the role of climate, hydrography and predation pressure on the seasonal long-term trends of each copepod. Each species' spring and summer biomass was then modeled by basin-specific smoothing functions of the different environmental variables, independent of whether long-term trend were significantly different between basins. The following model was used: with a b as the intercept at basin b and f b , g b , i b , and j b as basinspecific thin plate regression spline functions (maximum edf = 4) describing the effect of the locale climate index BSI, depth-specific temperature (T) and salinity (S), and predation (PI) respectively.
Prior to the analysis we tested for collinearity among the covariates using Pearson correlation coefficient and a variance inflation factor (VIF) of 3 as exclusion criterion [52]. In deciding which smoothing term to include in the final model, we applied a backward stepwise selection approach. Initially, we started with a full model that included a basin effect and basin-specific smoothers for all four covariates. As selection criterion we used the Akaike's Information Criterion (AIC) [53]. The underlying statistical assumptions were then tested through residual diagnostics of the optimal model. For the selected spring model of Acartia spp. residuals showed temporal autocorrelation of lag 2, hence we extended the GAM to a Generalized Additive Mixed Model (GAMM) by including a correlation structure. Such an extension allows a more complex stochastic model structure and implies that the single elements of the response variable are not independent anymore [51].We tested various structures in the full model and chose the one performing best based on the AIC [54]. We then applied again a model selection routine. All analyses were performed using the package 'mgcv' (with version R2.10) [51] within the 'R' statistical and programming environment [55].

Long-term trends
To understand how habitat heterogeneity affects the influence of climatic variability on zooplankton community structure and dynamics, we first examined the long-term development of our three indicator species Acartia spp., Temora longicornis, and Pseudocalanus acuspes in the Central Baltic Sea (CBS). We found T. longicornis and P. acuspes long-term trends to be spatially homogeneous, while long-term dynamics of Acartia spp. significantly differed between the three basins ( Table 1, Table S3).
The overall trend in Acartia spp. biomass was positive, particularly in spring in the Gdansk Deep (GD) where population sizes increased linearly with time (Fig. 2). This upward trend, however, leveled off at high biomass values at the end of the study period in the GD in summer as well as in both seasons in the Bornholm Basin (BB). In contrast, biomass in the Gotland Basin (GB) started to increase not until the late 1970s in spring with only a slight upwards trend in summer. T. longicornis biomass in spring increased as well, while in summer no significant trend was observed (Fig. 2, Table S3). P. acuspes biomass increased in all basins until the late 1970s, declined afterwards and increased again during the 2000s (Fig. 2).

External forcing
Significant external forcing variables for long-term trends of our zooplankton indicator species were strongly species-and seasonspecific, and partly basin-specific ( Table 2). The most parsimonious and best-fitting models for Acartia spp. spring and summer biomass explained 62% and 30% of the total variance, respectively. We found the climate index (BSI) to be the best predictor of annual variation in Acartia spp. biomass. For spring biomass the BSI effects were positive and linear, but of basin-specific strength (Fig. 3A). In contrast, the BSI effect on summer biomass was spatially consistent and non-linear with a stronger positive influence during periods with negative BSI anomalies (Fig. 3B). Temperature was a significant predictor of Acartia spp. biomass only in spring, but similar to the BSI of basin-specific strength. While the effect of temperature was weakly positive for the BB, it was highly significant and non-linear for the GD, with an overall negative trend.
Weaker predictors for Acartia spp. were salinity and predation, having uniform and generally negative linear effects in all basins. An exception represented the highly non-linear predation effect on spring biomass, being positive at the highest predation levels. As this is ecologically difficult to explain and rather points towards a spurious relationship, we conducted an alternative model selection excluding predation. This analysis resulted in significant autocorrelation, which we accounted for by including an exponential correlation structure. The selected alternative model contained only BSI as significant driver, but had less explanatory power than the model including predation (R 2 = 0.29, Fig. S3).
Our statistical analyses for T. longicornis revealed no basin effects ( Table 2). The most parsimonious and best-fitted spring model included, similarly to Acartia spp., a significant, albeit ecologically difficult to explain, predation effect (Fig. 3C). Hence we again applied an alternative model selection without predation, which resulted in the presented model, where most of the variance was explained by positive linear effects of temperature and BSI. Increasing salinity had additionally a strong, though negative influence. Differences between the models with and without predation were only minor, explaining 62% and 56% of the total variation, respectively. The variability in summer biomass, which lacked a clear trend, was only explained by the BSI. This driver showed a similar non-linear effect as with Acartia spp. but with limited explanatory power (explained only 5% of the fluctuations) and, thus, is not presented here.
Modeling long-term dynamics of P. acuspes revealed a particularly strong response to external drivers for summer biomass (the summer model explained 50% of the total variance) and no basinspecific relationships ( Table 2, Fig. 3D-E). The most important driver was predation pressure with a negative effect explaining 22% and 40% in spring and summer biomass variation, respectively. Temperature had a generally positive and slightly non-linear effect on spring biomass while salinity had an overall positive effect on the summer biomass. The BSI showed a clear negative linear effect on total biomass in spring, although this effect might be predominantly on the nauplii [36] which, like Acartia spp., inhabit the surface waters. In summer, a highly nonlinear BSI effect was observed with two difficult to explain maxima. Excluding this parameter would reduce the explanatory power of the model only by a few percent but would cause a strong impairment of the residual variance. We therefore kept the BSI in the final summer model.

Effect of spatial habitat heterogeneity
Our findings suggest vertical habitat heterogeneity strongly determines the effect of climate variability on the three most important zooplankton species with considerable effects on overall food web functioning. Acartia spp. inhabit the surface water layers in the central Baltic Sea and changes in this group were mainly affected by changes in atmospheric conditions (i.e., BSI). This finding agrees with previous studies reporting that changes in zooplankton phenology and biomass levels were particularly sensitive to variability in climate-driven physical factors during the beginning of the productive season [25,38,43]. The strong BSI effect on Acartia spp. compared to spring temperature suggests that other factors related to the BSI are of importance, such as vertical mixing intensity affecting the timing of the spring bloom and phytoplankton development [37,39]. Phytoplankton production is known to be an important driver for Acartia spp. [43], having a stronger effect on egg production than temperature [56,57]. Although Acartia spp. shows maximum egg production rates at Chl a concentrations of . 14-20 mg L 21 [57] the species can also establish dense populations with much lower food concentrations (,2 mg Chl a L 21 ) [58], suggesting that not only food quantity (with Chl a as proxy) but also food quality can be important for growth and reproduction (see review in [59]). Unfortunately, we lack sufficiently resolved long-term phytoplankton data to investigate the effects of food quantity and quality as well as seasonal phenology on CBS zooplankton dynamics.
Temora longicornis inhabits the mid-water habitat and changes in the biomass of this species were influenced by both atmospheric conditions and water temperature. The positive temperature effect on the long-term spring dynamics has been described before [25,38], whereas the relevance of primary production dynamics has only been shown recently [35]. Our data show that the longterm increase in spring biomass was, however, not translated to summer, as in Acartia spp., supporting the hypothesis that the increase reflects an earlier onset of population development rather than a numerical population response [35].
In contrast to Acartia spp. and T. longicornis, later developmental stages and adults of Pseudocalanus acuspes are directly exposed to fish predation in the deepwater habitat [31] explaining the statistical result of predation as the main driver of its long-term dynamics. The increase of the main predator sprat is due to increased temperature [60] and the decrease of its predator cod, the latter is partly due to lowered salinity and oxygen levels [61]. Hence, increased predation represents an indirect effect of climate [62]. P. acuspes showed also a strong positive response to increases in deepwater salinity in summer. This is in accordance with earlier studies claiming that this marine, glacial relict species prefers high salinities particularly for reproduction [25,33,63]. Younger stages (i.e., nauplii and copepodites C1-3) are highly abundant in spring, inhabiting the upper water layer [33]. Here they are affected mainly by temperature and other factors related to atmospheric forcing, such as phytoplankton dynamics [36].
In addition to the importance of vertical habitat heterogeneity for long-term dynamics of our study species, we were interested in the effect of spatial thermohaline variability on intra-specific Table 1. Summary of the F-ratio tests comparing species-and season-specific long-term trends modeled across basins or for each basin separately. population trends, which could have an effect on the community structure on a more regional scale. As this has not been adequately investigated in the realm of marine, pelagic systems, we compared long-term dynamics of populations of our study species that live in different basins of the CBS. We found the interaction of horizontal habitat heterogeneity with vertical habitat heterogeneity to be of importance. Although the CBS features a horizontal gradient in temperature and salinity, we did not observe basin-specific changes in biomass of T. longicornis or P. acuspes, which were strongly influenced by temperature and salinity respectively. A possible explanation is the temporal consistency of these gradients, as indicated by common long-term trends (Fig. S1). The lack of basin-specific responses indicates furthermore that differences in both hydrographic parameters seem not to be strong enough to cause basin-specific trajectories. For instance, instead of causing a much greater biomass decline of P. acuspes in the already less marine GB during the 1980s when salinity decreased, P. acuspes changed in the same way as in the BB, where salinity was always higher than in the GB (see Fig. S1). In contrast, Acartia spp. showed basin-specific long-term trends, coupled with spatial responses to changes in the BSI and temperature. The strong BSI effect on Acartia spp. biomass indicated that other factors such as wind mixing and consequently phytoplankton dynamics might be of importance. In spite of regional-scale atmospheric conditions wind-induced vertical circulation in the surface layers can show basin-specific patterns [64]. In line with this are the spatially differing long-term trends of different phytoplankton groups [39,65]. The resultant changes in diet composition could have affected the lipid dynamics and stoichiometric regulation of Acartia spp. and consequently the energy available for growth and reproduction [45,59].

Role of climate in shaping the ecological niche space
The importance of habitat heterogeneity for community structure and dynamics is tightly related to the concept of the ecological niche. In general, habitats with more structure allow the coexistence of species with differing physiological requirements and environmental preferences [17,66]. The range of a species' tolerance to the effects of multiple environmental factors was conceptualized by Hutchinson [67] as the n-dimensional hypervolume in environmental space that permits positive population growth, also termed as the fundamental niche (FN). This FN can also be understood as a fitness response-surface with an outer boundary defining the limit of population viability and inner contours representing increasing fitness [68]. Often, in the real world, not all combinations of environmental conditions, which could be favorable, are realized at a particular time or occur within a region. Consequently, only a subset of the FN exists which has been termed the potential niche (PN) by Jackson and Overpeck [69]. The PN may be substantially smaller than the FN [70]. The realized niche (RN) in this context represents the part of the PN that is constraint by biotic interactions (e.g., competition, predation) [67] or enlarged under high species dispersal [71,72]. Hence, the RN can be envisioned as the populations' response-surface in terms of, e.g., density or biomass [73]. When climate changes, the PN of a species can change in shape or size and similarly its RN due to additional changes in biotic interactions. Some species may benefit from these changes and increase in abundance, other species may persist but reduce in numbers, migrate or undergo local extinctions [74,75]. Consequently, changes in community structure at a given locality depend not only on the magnitude of environmental changes but also on the level of species-specific niche differentiation.
In this respect, the increase in population size of Acartia spp. and T. longicornis can be interpreted as an improvement in the thermal space of their PN since both species showed a positive response to increasing spring temperatures. Both species further show highest egg production rates at similar temperatures, i.e., 13-18uC for  Acartia spp. and 17uC for T. longicornis, and high mortalities at 24uC [32,57], suggesting a great overlap in the optimal thermal space of their FN. Consequently, the PN of these two species in the upper 40 m of the CBS was likely sub-optimal in spring since water temperatures ranged generally only between 1 and 8uC (Fig. S1). A further warming would likely lead to increased productivity of these two species. This could be the case particularly for T. longicornis as this species responds stronger to temperature changes than Acartia spp. [76]. The role of biotic interactions such as predation becomes particularly apparent in the long-term dynamics of P. acuspes. Little is known about the response of its physiological processes to synergistic effects of environmental factors and, consequently, the size or shape of its FN or PN. But the fact that we identified predation as the most important factor for variation in the populations' biomass suggests that mainly the RN space changed over time. For the marine copepod Calanus finmarchicus, Helaouët & Beaugrand [77] observed a strong coupling of the RN and the FN with rather similar niche widths. This congruency was explained by limited migration due to hydrodynamical barriers while biotic interactions were not discussed. For Baltic Sea P. acuspes, however, the RN is likely to be smaller owing to a strong top-down control through planktivorous fish. In general, the level of predation pressure can be reduced by behavioral adaptation. For instance, optical sampling has recently revealed predator avoidance behavior of P.acuspes females in the Baltic, which forces a high proportion of the population into sub-optimal habitats [78]. These habitats below the permanent halocline have low dissolved oxygen concentrations and are avoided in the absence of those predators. Such escapement behavior therefore only extends the outer boundary of the RN. In other words, the space in which the population is viable increases but not the space for optimum fitness [68].

Conclusions
Our study findings suggest that habitat heterogeneity can modulate the response of zooplankton communities to climate changes and the extent in which the community structure can be altered. In the Baltic Sea, habitat heterogeneity has promoted a mixture of copepod species to exist in the central basins that have different fundamental niches and that display different sensitivities to climate-driven forcing. Consequently, the long-term development of these species differed greatly under changing environmental conditions, which in turn influences the community structure and dynamics. Individual habitats may also respond differently to climatic variability, as seen for the surface layers in the different basins, leading to differences in population trends even within the same species. We stress therefore the importance of understanding how anticipated climate change will affect ecological niches and habitats in order to project spatio-temporal changes in species abundance and distribution. Figure S1 Long-term trends of the hydrography within the Central Baltic Sea (CBS). The main trend of the spring and summer temperature and salinity in different water layers (see Table S2) are given for each basin (red = Bornholm Basin BB, blue = Gdansk Deep GD, green = Gotland Basin GB) or the entire CBS region if trends and overall values were basin-unspecific (black line). Colored points represent the observed values in each basin, while the continuous lines indicate the predicted trends from Generalized Additive Models (GAM). Grey shaded areas indicate the associated standard errors.