Long-term changes in spatial overlap between interacting cod and flounder in the Baltic Sea

The strength of interspecific competition and predator–prey interactions depends on the area of co-occurrence of the interacting species. Therefore, it is necessary to quantify the changes in the spatial overlap of trophically connected species to understand the outcomes of species interactions. In the Baltic Sea, the interplay between cod and flounder has previously been neglected. In this study, we use four decades of data on cod and flounder distributions covering the southern and central Baltic Sea to: (1) model and map the changes in the distributions of the two species using generalized additive models; (2) quantify the temporal changes in the potential competitive and predator–prey interactions between them using spatial overlap indices; (3) relate these changes in overlap to the known dynamics of the different cod and flounder populations in the Baltic Sea. Competition overlap has continuously increased for cod, from the beginning of the time-series. This is a possible cause of the observed decline in feeding levels and body condition of small and intermediate sized cod. Flounder overlap with large cod instead has decreased substantially, suggesting a predation release of flounder, potentially triggering its increase in abundance and distribution range observed in the last decades.


Introduction
A key aim of ecology is to understand how species interact with each other and how their population dynamics are shaped because of these interactions. Changes in interspecific interactions through predation and resource competition can affect for example abundances, spatial distributions, mortality rates, growth, fecundity, and age/size at first maturity (Sparholt, 1994;McPeek & Peckarsky, 1998;Klanderud, 2005;Tylianakis et al., 2008;Hillyer & Silman, 2010). In aquatic systems, these effects may in turn produce massive reorganizations of the food webs by means of trophic cascades resulting in regime shifts (Gårdmark et al., 2015;Rudman et al., 2016;Walsh et al., 2016), with serious implications for species conservation and the management of exploited populations. In particular, quantifying the spatial dynamics of species distributions has proven to be a central aspect for understanding the outcomes of species interactions, including how spatial heterogeneity affects functional responses, or how the results of small-scale experimental studies on trophic interactions can be extrapolated to heterogeneous natural systems (Chesson, 1998;Englund & Leonardsson, 2008).
The spatial distribution of species is influenced by environmental factors, as well as by interactions with other species and anthropogenic activities (Hutchinson, 1957;Wisz et al., 2013). Moreover, the spatial and temporal patterns in species distributions as well as their overlap determine the strength of trophic interactions such as predation and competition (Chesson, 1998;Garrison, 2000;Bergström et al., 2006;Englund & Leonardsson, 2008). Species distribution models relate species distribution data (occurrence or abundance) to environmental explanatory variables (Elith & Leathwick, 2009) and have been widely used to map spatial distributions of species and populations and to predict changes in response to an altered environment (Austin, 2002;Elith & Leathwick, 2009;Bergström et al., 2013;Pellissier et al., 2013). In aquatic systems, for example, the distribution maps obtained from these models have been used to estimate overlap between species, populations and life stages, which can in turn be used to infer interaction strength (Hunsicker et al., 2013).
In the brackish and species-poor Baltic Sea, cod (Gadus morhua, Gadidae) and the flounder species complex (Platichthys flesus and Platichthys solemdali, Pleuronectidae, hereafter simply referred to as flounder, Momigliano et al., 2018) are the main demersal fish species, in terms of abundance and commercial importance (Casini et al., 2008;Florin & Höglund, 2008;Lindegren et al., 2009). The Eastern Baltic cod (hereafter simply referred to as cod) and flounder populations have undergone major changes in their abundance, spatial distribution and growth in the last century (Eero et al., 2007;Casini et al., 2012;Bartolino et al., 2017;Orio et al., 2017aOrio et al., , 2019. Trophic interactions between cod and flounder in the Baltic have been suggested to play a role in these changes (Persson, 1981;Orio et al., 2017aOrio et al., , 2019, since large cod individuals predate on flounder (Almqvist et al., 2010;ICES, 2016) and cod and flounder may compete for benthic prey (Arntz & Finger, 1981;Gjøsaeter, 1988;Haase, 2018). However, no formal analyses of spatial overlap have been performed so far that could shed light on the potential interactions between cod and flounder in support to these hypotheses.
Since the strength of species interactions depends, besides species densities, also on the area of their cooccurrence (Chesson, 1998;Bergström et al., 2006), it is necessary to quantify the changes in the spatial overlap of trophically connected species (Hunsicker et al., 2013). Therefore, in this study, we first model and map the changes in the distributions of different size classes of cod and flounder in the Baltic Sea. We then use different overlap indices to illustrate the changes in the potential competition and predatorprey interactions between these species both on a spatial and a temporal scale. Finally, we relate these changes in overlap to the known dynamics of the different cod and flounder populations in the Baltic Sea and provide a practical example on how our results could be used in further analyses.

Materials and methods
Trawl survey data Standardized catch per unit of effort (CPUE) in weight (kg h -1 ) of cod and flounder were obtained from a database containing data collected during the Baltic International Trawl Survey (BITS; ICES 2014) and historical Swedish and Latvian trawl surveys carried out between 1979 and 2016 in ICES Subdivisions (SDs) 22-29 ( Fig. 1, Online Resource Appendix 1, Fig. A1, Table A1). The BITS survey has cod and flatfish as target species (ICES 2014) and the data from this survey are routinely used in the stock assessments and to provide fisheries advice on these species in the Baltic Sea (ICES 2019).
For information about the standardization and the survey database, see Orio et al. (2017a). As in Orio et al. (2019) we included in the analyses all the hauls performed in SDs 24 to 28, excluding the Gulf of Riga (SD 28-1), in quarters 1 and 4 with depths between 20 and 120 m because of good spatial and temporal coverage. We excluded from the analyses all the hauls performed in SD 27 north of 58 degrees latitude and in SD 24 west of 12.5 degrees longitude, as these areas were not consistently covered by the bottom trawl surveys during the period analyzed.
We aimed at studying spatiotemporal changes in the potential interactions between cod and flounder in terms of both competition and predation. As in Orio et al. (2017aOrio et al. ( , 2019 we assumed that the spatiotemporal changes in the total flounder CPUEs would reliably represent the trends of the larger individuals of the population (C 20 cm length). Here we assume that both competition and predation act in the same way for the two flounder species since they are identical in shape and their diet comprises mostly benthic food (Momigliano et al., 2018). To investigate the competition between cod and flounder we focused only on cod between 15 and 35 cm because at this size zoobenthos constitutes a major proportion of their diet in the Baltic Sea (ICES, 2016), as is the case also for large flounder (Florin, 2005;Haase, 2018). To investigate the predator-prey interactions we considered only cod C 55 cm since from around this length cod predate on flounder bigger than * 20 cm in the Baltic Sea (Almqvist et al., 2010, ICES, 2016).

Generalized additive model (GAM)
The CPUEs of cod and flounder were transformed in presence/absence data and modeled using generalized additive models (GAM) using a binomial error distribution with a logit link function (Wood, 2006). Three models were built to explain the distributions of cod 15-35 cm, cod C 55 cm , and flounder.
The full binomial models for presence/absence were formulated as follows: where b is an overall intercept different for each quarter, s is an isotropic smoothing function (thinplate regression spline; Wood, 2003), te is a tensor product smoothing function used for representing interaction terms, f i are natural cubic splines, and e is an error term. The interactions were introduced to take into account the changes in the spatiotemporal distribution of the species over the long time period analyzed. Data on bottom salinity, temperature, and oxygen were monthly averages of the selected years extracted from the hydrodynamic Kiel Baltic Sea Ice-Ocean Model (BSIOM; Lehmann & Hinrichsen, 2000, Lehmann et al. 2002. Bottom temperature, salinity, and oxygen were included in the analyses to increase the predictive power of our models, since these variables have been shown to be important in explaining the distribution of both gadoids and flatfishes (Begg & Marteinsdottir, 2002;Hedger et al., 2004;Able et al., 2005;Hinrichsen et al., 2009).
Model selection was done through a backward stepwise elimination process based on statistical significance (Wood, 2006). From the full model, the nonsignificant predictor with the lowest significance level was excluded at each step and the model run again. This procedure was repeated until all the predictors were significant (final model; Orio et al., 2017aOrio et al., , 2019. To obtain ecologically significant models and to avoid overfitting, we set a limit to the maximum degrees of freedom (number of knots, k) allowed to the smoothing functions of the variables latitude, longitude, depth, temperature, salinity and oxygen (k = 4) and of the interaction between latitude and longitude (k = 20).
Multi-collinearity in the models was checked using the variance inflation factors (VIF) calculated with the usdm library of R (Naimi et al., 2014). Spatial autocorrelation was checked using spline correlograms produced with the ncf library of R (Bjornstad, 2018), while temporal autocorrelation has been assessed using the acf function in R (R Core Team 2017).

Mapping probability of co-occurrence
The final models were used to predict probability of occurrence of cod 15-35 cm, cod C 55 cm and flounder over a regular grid of 0.02°9 0.01°(corresponding to approx. 1 9 1 km) for quarter 1 and quarter 4. The same areas and depths removed from the data used for modeling were removed from the prediction grid. The environmental data associated with the grids were the averages of the first 3 months of each year for quarter 1 (January-March) and the last three for quarter 4 (October-December). The predicted probability of occurrence for flounder was then multiplied by the predicted probability of occurrence for both cod 15-35 cm and cod C 55 cm, respectively, to obtain the predictions of the probability of co-occurrence for both competition and predator-prey interaction.

Reconstructing changes in spatial overlap
In order to quantify potential changes in trophic interactions, we calculated an index of spatial overlap that takes into account both the cod and the flounder perspectives. Following the approach of Neuenfeldt & Beyer (2006), from the annual spatial predictions of each final model we calculated the percentage of area where species A occurs with a probability of occurrence C 0.75 in which also species B occurs with a probability of occurrence C 0.75, and vice versa. This approach estimates the relative encounter space between the two species. For predator-prey interactions, this approach tackles the two operational ways to formulate predator-prey overlap (i.e., from the predator and prey perspectives, quantifying the prey availability for the predator and the predation risk for the prey, respectively). In this way, we obtained timeseries of percentage of area in which there is potential predation or competition between the two species from both the cod and flounder perspectives. Beside co-occurrence, species densities are also important in the realization of an interaction, and therefore we explored the potential differences in the overlap that could be obtained by using abundance data. To this end, we ran Delta-GAMs as in Orio et al. (2017aOrio et al. ( , 2019 and modified the overlap indices as described in the Online Resource Appendix 1. All the analyses were performed using the R software and the mgcv and ggplot2 libraries of R (Wood, 2011;Wickham, 2016;R Core Team, 2017).

Results
The analyses included a total of 9969 hauls for cod and 9764 hauls for flounder. The final binomial models explained 20.1-31.6% of the deviance depending on species and size class ( Table 1). The high degree of unexplained variation, which is common in large-scale fish models (Grüss et al., 2014;Parra et al., 2016), probably originates from noise introduced by the fish sampling method and the use of modeled environmental data, in combination with stochasticity in fish distributions ). Thus, we believe that the models still captured the main general patterns in the data. Since the VIF values were B 3 (Zuur et al. 2009), multi-collinearity was considered negligible. Residuals of the models were visually inspected and revealed, in some instances, slight departures from the model assumptions, but the general quality of the residuals was deemed satisfactory (Online Resource Appendix 1, Figs. A2-A4). No spatial autocorrelation of the residuals of the models based on 95% pointwise bootstrap confidence intervals was found after plotting spline correlograms (Online Resource Appendix 1, Figs. A5-A7). Small but significant temporal autocorrelation was detected but, given the nature of the dataset used and the predictive scope of the models, it was considered acceptable.
All the interactions were retained in the final models. Regarding the variables not in interactions, only oxygen, bottom temperature and quarter were retained in all models. The partial effects of the final models are presented in the Online Resource Appendix 1, Figs. A8-A10.
The predictions of the spatial distribution of the probability of occurrence of cod 15-35 cm, cod C 55 cm and flounder over time are presented in Online Resource Appendix 1 Figs. A11-A13.
The predicted probability of co-occurrence (Figs. 2 and 3) were higher and more widespread in the first quarter of the year compared to the fourth quarter, for both competition and predator-prey interaction. However, for competition (i.e., between cod 15-35 cm and flounder; Fig. 2) the probability of co-occurrence increased throughout the period investigated, while for predator-prey interaction (i.e., between cod C 55 cm and flounder; Fig. 3) it generally decreased.
The time-series of potential competition from the cod perspective (percentage of area where cod 15-35 cm is present and in which also flounder occurs; Fig. 4 top panels) showed values lower than 50% at the beginning of the time-series in all SDs n numbers of stations used in the models, edf effective degrees of freedom, Dev% explained deviance Fig. 2 Predictions of the probability of co-occurrence of cod 15-35 cm and flounder (competition) in the first and fourth quarter of selected years across the study period except for SD 24 and then a rapid increase from the mid-1980s up to * 80-100% depending on the quarter. In general, the area of potential competition for cod was larger in the first quarter. From the flounder perspective, instead, the time-series of potential competition (percentage of area where flounder is present and in which also cod 15-35 cm occur; Fig. 4 bottom panels) showed values close to 100% in the first years, then a decrease with a minimum around the mid-1990s (down to 10-20% in the northern SDs) and a subsequent re-increase. The amplitude of this U-shaped trend was larger in the northern SDs. In the southern SDs, moreover, the area of potential competition for flounder reached in the most recent years values as high as at the beginning of the timeseries. In general, the area of potential competition for flounder was larger in the fourth quarter. The time-series of potential predation from the cod perspective (percentage of area where cod C 55 cm is present and in which also flounder occurs, representing prey availability for large cod; Fig. 5 top panels), showed, as for competition, values lower than 50% at the beginning of the time-series in all SDs except in SD 24 and then a rapid increase from the mid-1980s up to * 80-100% depending on the quarter. However, in the northern SDs cod C 55 cm disappeared from the beginning of the 1990s as shown by the broken lines in Fig. 5a. In general, the prey availability (i.e., the availability of flounder as prey) for large cod was higher in the first quarter. From the flounder perspective, instead, the time-series of potential predation (percentage of area where flounder is present and in which also cod C 55 cm occurs, indicating predation risk for flounder; Fig. 5 bottom panels) showed values close to 100% in the first years and then a drop until recent years with only a small peak around the late-2000s in the southern SDs. In the northern SDs, the predation risk of flounder had been close to 0 from the beginning of the 1990s. The trends in predation risk were very similar in the two quarters.
The time-series of area of overlap obtained from the abundance estimates are relatively similar to those based on presence/absence data (Online Resource Appendix 1, Figs. A14-A15), with the exception of Fig. 3 Predictions of the probability of co-occurrence of cod C 55 cm and flounder (predator-prey interactions) in the first and fourth quarter of selected years across the study period the percentage of area with a potential competition from the flounder perspective in SD 25, which shows an increase from the beginning of the time-series instead of a U-shaped trend (Online Resource Appendix 1, Fig. A14). Therefore, we are confident that the trends presented using presence/absence data reliably picture the changes in the potential ecological interactions between cod and flounder.

Discussion
Our analyses reveal pronounced changes in the spatial overlap and in the potential competitive and predatorprey interactions between cod and flounder in the Baltic Sea in the last four decades.
In general, the maps show clear seasonal differences with higher and more widespread probabilities of co-occurrence in the first quarter of the year compared to the fourth quarter. This result is likely due to the fact that more flounder is found in offshore areas in the first quarter compared to the fourth, as shown in Fig. A13. Cod, instead, show a less distinct seasonal difference in distribution (Figs. A11 and A12). The more widespread distribution of flounder during the first quarter of the year can be explained by the fact that flounder's spawning season starts at the beginning of spring (especially in the case of the pelagic spawning flounder P. flesus; Orio et al., 2017b). Therefore, the trend shown in the first quarter could be linked to the spawning migration of flounder toward deeper areas.
The maps show high probabilities of co-occurrence of cod 15-35 cm and flounder throughout the period investigated and an increase toward the more recent years. On the contrary, the probabilities of co-occurrence of cod C 55 cm and flounder declined during the period investigated with a particularly intense decline in the northern areas. However, the time-series of percentage of area with a potential competition and predator-prey interaction between these species

Potential mutual competition between cod and flounder
The potential competition for cod (i.e., the area with cod where also flounder occurs; Fig. 4 top panels) was low at the beginning of the time-series and increased steeply from the mid-1980s, while the potential competition for flounder (i.e., the area with flounder where also cod occurs; Fig. 4 bottom panels) shows a U-shaped trend with high values at the beginning and end of the time-series and lower values in the mid-1990s. These changes in potential competition between cod and flounder are in line with their known population dynamics. In the beginning of the 1980s, cod in the Baltic experienced a massive increase in abundance and it was widely distributed in areas where it is usually not present (Casini et al., 2012;Casini, 2013;Orio et al., 2019). On the contrary, flounder abundance was low and its distribution was concentrated in a smaller portion of the Baltic . These dynamics can explain the low potential competition from the cod perspective and the high potential competition from the flounder perspective at the beginning of the time-series. The increase in potential competition for flounder at the end of the time-series can instead be explained by the changes in the depth distributions of both cod and flounder toward more similar depths, likely due to increased deepwater hypoxia . In the case of cod, its collapse, the contraction of its geographical and depth distributions and the following increase in abundance and distribution of flounder  are the underlying processes that caused the steep increase in potential competition with flounder in the end of the 1980s revealed by our study. During the same period the body condition of cod started decreasing and stabilized at very low levels around the beginning of  (Casini et al., 2016). This drop in cod body condition has been related in literature to the lack of pelagic prey and the increase of hypoxic area (Casini et al., 2016). Furthermore, recent analyses revealed a decreased feeding level for cod \ 30 cm from the beginning of the 2000s, mainly due to a decrease in benthos feeding (ICES, 2017a, Neuenfeldt et al., 2020. This could be explained by a decrease in benthic prey availability as a consequence of the massive increase in hypoxic areas that occurred after the mid-1990s (Casini et al., 2016). Here we argue that the increase in competition for benthic food with flounder could have contributed to the low feeding level and the drop in body condition of cod. Competition is probably directed to especially the benthic isopod Saduria entomon, that is an important prey for both cod and flounder in the offshore Baltic Sea (Haase, 2018). To test this hypothesis and to provide a concrete example of how the overlap indices we produced can be used in further analyses, we performed a correlation analysis between the potential competition from the cod perspective (as estimated by our study) and the body condition of cod 20-30 cm in the fourth quarter in SDs 25-28 (updated from Casini et al., 2016). The analysis shows a temporal negative correlation (r = -0.61; Fig. 6) supporting this hypothesis.
The hypothesis of high competition for benthic food between small cod and flatfishes was previously also suggested by Persson (1981) to explain the low abundances of cod at the beginning of the twentieth century. From Fig. 6, it is also possible to note that considering only the years from 1985 the relationship between area of overlap and cod condition becomes stronger (r = -0.75). In the mid-1980s, a major regime shift occurred in the Baltic with a massive reorganization of the Baltic food-web structure and functioning Casini et al., 2009). Such a regime shift could explain the two periods identified in Fig. 6. In the first half of the 1980s, the area of overlap between cod and flounder was not an important factor determining cod condition, possibly due to the low biomass of flounder, while, from the mid-1980s, the increase in area of overlap corresponds to a decrease in body condition. In the last 20 years, however, body condition of small cod kept decreasing fast although the area of overlap declined at a lower rate, which could be a sign that other factors are strengthening the effect of competition between cod and flounder, as for example the increase in anoxic areas causing a further decrease in the benthic fauna, or lack of pelagic prey (Casini et al., 2016).

Cod prey availability and flounder predation risk
The percentage of area with flounder prey available for large cod (i.e., the area with large cod where also flounder occurs) increased from the mid-1980s. This could partially explain the results of Neuenfeldt et al. (2020) and ICES (2017a) who found that from the mid-1990s the feeding levels of cod C 55 cm have increased and are the highest recorded from the 1960s. However, the effects of the high feeding levels of large cod were not reflected in high body condition (Casini et al., 2016) potentially due to a combination of different processes, other than food availability, affecting cod condition (Orio, 2019). For example, the increased intensity and prevalence of parasite infection in cod (ICES, 2017b), especially the large ones, could be a process that may affect energy conversion in cod (Horbowy et al., 2016). Moreover, it has been proven that hypoxia can alter the metabolism and growth of organisms (Chabot & Dutil, 1999;Diaz & Rosenberg, 2011;Levin, 2018). Therefore, cod living in areas with low oxygen concentrations could experience physiological stress affecting its condition even if its feeding level is high (Orio, 2019). Also, the decrease in benthic preys due to the increase in Fig. 6 Correlation between body condition of cod 20-30 cm (updated from Casini et al., 2016) and percentage of area in which there is potential competition between cod 15-35 cm and flounder from the cod perspective in the fourth quarter in SDs 25-28. The line connects subsequent years hypoxic areas could have forced cod to predate relatively more often on pelagic species, which requires higher energy, thus potentially contributing to explain cod low body condition. It is also important to consider that the abundance of cod C 55 cm has decreased drastically from the mid-1990s, and in the most recent year they are very seldom caught (Orio, 2019). Therefore, a definite conclusion on what is affecting the feeding level of such a poorly represented length class of cod cannot be drawn.
From the flounder perspective, instead, the predation risk from cod has decreased steadily in all areas, especially in the north (SDs 27 and 28), following the dynamics of the cod stock, which during the same period collapsed and contracted to the southern areas of the Baltic (Bartolino et al., 2017;Orio et al., 2017aOrio et al., , 2019. Moreover, cod maximum length also dropped at the same time (Orio et al., 2017a), further decreasing the amount of large cod capable of feeding on flounder. These changes in the cod stock strongly suggest a predation release of flounder that favored the observed increase in its abundance and extent of distribution (Orio et al., 2017a, as also happened for small pelagic species both in the Baltic and in other areas where cod stocks collapsed (Frank et al., 2005;Casini et al., 2009).

Interactions between gadoids and flatfishes
Diet overlap and potential competition between gadoids and flatfishes have been studied across the Atlantic Ocean. Similar diets between gadoids and flatfishes have been recorded for example both on the Scottish coast (Gibson & Ezzi, 1987) in the NE Atlantic, and at Georges Bank (Link et al., 2005 and reference therein) in the NW Atlantic. In these studies, competition has been advocated as a hypothesis to explain trends in abundance of flatfishes and other fish species including gadoids. However, in the case of the Scottish coast (Gibson & Ezzi, 1987), food resource competition between flatfishes and gadoids was only speculated since no information on resource abundance was available. In the Baltic Sea case, there are indications that currently the benthic resources could be limited due to the large extent of hypoxic areas (Karlson et al., 2002;Casini et al., 2016;Karlson et al., 2019). In fact, the spreading hypoxia has caused a lack of macrofauna on large areas of the sea bottom during the past two decades, with repercussions on all trophic levels of the Baltic Sea ecosystem (Karlson et al., 2002;Villnäs et al., 2013). Specifically, according to Eero et al. (2012) the lack of benthos could have been one of the causes triggering density-dependent effects in cod, such as increased cannibalism and decreased growth, while for Neuenfeldt et al. (2020) it could be linked to the decreased feeding level for cod \ 30 cm.
Gadoid predation on flatfishes has also been recorded in other areas of the North Atlantic, such as the North Sea, the Scottish coast and the Bering Sea (Gislason & Helgason, 1985, and references therein, Bailey, 1994and references therein, Ellis & Gibson, 1995, Van der Veer et al., 1997. In the Bering Sea, flatfish predation on gadoids has also been described and was speculated as having a potential effect on the recruitment of gadoid populations (Hunsicker et al., 2013). In the Baltic, however, there is no indication that flounder predate on small cod and therefore affect its recruitment success.

Conclusions
In the Baltic Sea, studies of fish species interactions in the offshore have focused mainly on the interplay between cod, sprat, and herring. Our study provides for the first time insights into the potential interspecific interactions between cod and flounder in the Baltic Sea in the last four decades, which could contribute to explain some of the observed changes in their population dynamics. Competition overlap has increased for cod from the beginning of the timeseries, constituting a possible cause of the low feeding levels and body condition of small and intermediate sized cod, in addition to the increase in anoxic areas and decline in pelagic prey (Casini et al., 2016). Predation risk instead has decreased substantially for flounder, potentially triggering its increase in abundance and extent of distribution in the Baltic. The investigation of species interactions by means of species distribution models and overlap indices, as shown in this study, offers a valuable set of complementary data to the more common information on species interaction based on dietary analyses and can be used to estimate population-level effects of interactions in heterogeneously distributed species (Englund & Leonardsson, 2008). Our results are also highly relevant for ecosystem-based fisheries management as they can be implemented, for example, in multispecies models to inform on the effects of spatial heterogeneity on the functional responses of interacting species.