Substantial changes in the depth distributions of benthic invertebrates in the eastern Kattegat since the 1880s

Abstract Bottom trawling and eutrophication are well known for their impacts on the marine benthic environment in the last decades. Evaluating the effects of these pressures is often restricted to contemporary benthic data, limiting the potential to observe change from an earlier (preimpact) state. In this study, we compared benthic species records from 1884 to 1886 by CGJ Petersen with recent data to investigate how benthic invertebrate species in the eastern Kattegat have changed since preimpact time. The study shows that species turnover between old and recent times was high, ca. 50%, and the species richness in the investigation area was either unchanged or higher in recent times, suggesting no net loss of species. Elements of metacommunity structure analysis of datasets from the 1880s, 1990s, and 2000s revealed a clear change in the depth distribution structure since the 1880s. The system changed from a Quasi‐nested/Random pattern unrelated to depth in the 1880s with many species depth ranges over a major part of the studied depth interval, to a Clementsian pattern in recent times strongly positively correlated with depth. Around 30% of the 117 species recorded both in old and in recent times, including most trawling‐sensitive species, that is large, semiemergent species, showed a decrease in maximal depth of occurrence from the deeper zone fished today to the shallower unfished zone, with on average 20 m. Concurrently, the species category remaining in the fished zone was dominated by species less sensitive to bottom trawling like infauna polychaetes and small‐sized Peracarida crustaceans, most likely with short longevity. The depth interval and magnitude of the changes in depth distribution and the changes in species composition indicate impacts from bottom trawling rather than eutrophication. Furthermore, the high similarity of results from the recent datasets 10 years apart suggests chronic impact keeping the system in an altered state.


| INTRODUC TI ON
Almost no areas globally are now unaffected by humans, and many areas, including the Scandinavian seas, are regarded highly impacted by several pressures (Halpern et al., 2008(Halpern et al., , 2015. Therefore, a frequent problem when evaluating human environmental impact is the paucity of good contemporary reference data. Two major recent threats to the marine fauna and flora in Scandinavian coastal marine environments are eutrophication and fishery impacts due to bottom trawling, recognized in the recent 3-4 decades. Here we take advantage of a temporal reference, the benthic invertebrate records from the surveys with the canon ship "Hauch," performed in the 1880s in the Kattegat by CGJ Petersen (1893). This is the earliest comprehensive survey dataset on benthic invertebrate fauna available from the Scandinavian Sea area, from a time before major impact of the mentioned threats. The surveys took place in the beginning of the industrial expansion in the Scandinavian countries, before the rational agriculture and the application of water closets, and before intensive use of bottom trawls, which required strong motor-driven vessels. Effective trawling therefore increased sharply first in the 1930s with the introduction of diesel engines (Lindeboom & de Groot, 1998).
Several previous studies have used benthic historic data to set reference conditions in the Scandinavian Sea area. Stephenson, Williams, and Cook (1972) reanalyzed Petersen data from the Danish waters taken some 25 years later than the data from Petersen (1893), and the binary approach mainly confirmed the broad categorization into "Petersen communities." Revisits of historic sites in the 1970s and 1980s focusing on biomasses reported some support for organic enrichment of the area. Cederwall and Elmgren (1980) visited sites in the Baltic Sea, Pearson, Josefson, and Rosenberg (1985) sites in the Danish sector of the Kattegat established in 1911-1912(Petersen, 1913, and Rosenberg, Gray, Josefson, and Pearson (1987) sites in the Skagerrak and Oslofjord from 1914 (Petersen, 1915). As the Petersen investigations in 1911Petersen investigations in -1914 had the purpose to determine the amount of fish food, the taxonomic resolution was much lower than in the study in the 1880s. For instance, polychaete species were often lumped into one category, "vermes." To our knowledge, our study is the first one using historical data to examine changes in spatial distributions of species in the Kattegat in relation to different environmental threats.

Effects of eutrophication have been obvious in the
Scandinavian sea areas during the second half of the twentieth century, an impact now showing signs of reversal after peaking in the 1980s (Andersen et al., 2017;Carstensen, Conley, Andersen, & AErtebjerg, 2006;Riemann et al., 2016). Concurrently, there are reports of decreased Secchi depths in recent times (e.g., Middelboe & Sand-Jensen, 2000), which, at least partly, may be related to eutrophication. This together with reports on deeper distribution of eelgrass (Zostera marina) in the late 1800s (Petersen, 1893) than in the second half of the 1900s (e.g., Duarte, 1991) indicates a narrowing of the primary productive euphotic zone, with potential effects on benthic micro-and macroalgal consumers like grazing snails and some filtering bivalves. Another effect of eutrophication, predominant in some stagnant deeper areas, is the promotion of hypoxia/anoxia, causing mortality of benthic invertebrates (Diaz & Rosenberg, 2008).
Several fishery activities have impacts on marine ecosystems (Dayton, Thrush, Agardy, & Hofman, 1995). While impact due to overfishing, particularly removal of large top predators, has occurred for a very long time (Jackson et al., 2001), bottom trawling is a more recent phenomenon and poses a particular threat to benthic habitats (Thrush & Dayton, 2002). In addition to direct injury of benthic organisms, benthic trawling may homogenize the benthic environment and thereby reduce the habitat complexity at several scales, with adverse effects on biodiversity (Hewitt, Thrush, Halliday, & Duffy, 2005;Thrush & Dayton, 2002;Thrush et al., 1998;Watling & Norse, 1998). Several recent studies from the northeastern Atlantic area report adverse effects of trawling on the benthic fauna such as reduced biodiversity, biomass, and production (Hiddink et al., 2006;Kaiser & Spencer, 1996;Sköld et al., 2018;Van Denderen, Hintzen, Rijnsdorp, Ruardij, & van Kooten, 2014). Trawling may also affect species trait composition by causing a shift from communities dominated by relatively sessile, emergent, high biomass species to communities dominated by small-bodied infauna (Kaiser, Ramsay, Richardson, Spence, & Brand, 2000). In the area of the present study, effects of trawling are likely depth dependent because intense trawling activities nowadays mainly take place at water depths greater than ca. 30 m (ICES, 2016;Pommer, Olesen, & Hansen, 2016).
Because potential effects of the mentioned threats are depth dependent, our main hypothesis is that the depth distribution structure of the benthic species has changed since preimpact time in the 1880s. We apply the elements of metacommunity structure (EMS) approach (Leibold & Mikkelson, 2002), an unbiased method where the data determine the pattern. The method has proven effective to identify the distribution pattern along the main gradients in several realms (Heino, Soininen, Alahuhta, Lappalainen, & Virtanen, 2017;Josefson, 2016;Meynard et al., 2013;Presley, Higgins, & Willig, 2010;Valanko, Heino, Westerbom, Viitasalo, & Norkko, 2015). We expected trawling-sensitive species (Figure 1) to disappear or occur only in nonfished areas away from the trawled zone if trawling impact was important, leaving richness in trawled areas dominated by trawling-insensitive species.

| Sampling
We compare the water depth distribution structure and species richness in the 1880s dataset with datasets from the 1990s and the 2000s in the same subarea (i.e., the Swedish sector of the Kattegat) at sites with similar bottom type and in the same depth interval, that is 19-80 m (Figure 2, Supporting Information Table S1).
The 1880s: The 1880 dataset is qualitative (binary) with the aim to describe the spatial distribution of the common species in the area (Petersen, 1893). Consequently, the 1880 data from the whole Kattegat have a high spatial (a total of >500 sites) and a high taxonomic resolution (>570 species).
We selected 57 sites in the depth interval of 19-80 m from the Swedish sector of Kattegat for study, and the bottom substrate was at least partly fine sediments like sand, silt, or mud (Figure 2, Supporting Information Table S1). We determined positions of sampling sites ourselves from bearings and distance from fixed points like lighthouses and water depths given in fathoms (Petersen, 1893), seemingly in fair agreement with recent bathymetrical maps. Sampling in the 1880s used both trawls and dredges, mostly simultaneously at the same site. There is no information on mesh sizes of gears or of eventual sieves for extraction, but figures given on body sizes of a few millimeters of many mollusk species and records of many small-bodied crustaceans like most cumaceans indicate a high degree of detail in the examination of the samples (Petersen, 1893). The surveys took place in the four sequential years 1883-1886. Extraction of fauna from these samples was performed using 1-mm sieve and otherwise followed standard procedures in the Baltic Sea area (HELCOM, 1988).

| Data treatment and approach
Taxonomic identities were compared and adjusted using the species list WorMS http://www.marinespecies.org/.
We applied community analyses on the full datasets from each decade. Additionally, in order to evaluate whether the depth distribution of the same individual species had changed over the >100 year period, analyses were run on the subset of species recorded both in old and in recent times (Supporting Information Tables S2 and S3).

| Elements of metacommunity structure analysis
The elements of metacommunity structure (EMS) approach uses a stepwise procedure and can simultaneously test for multiple idealized patterns across a set of sites (Leibold & Mikkelson, 2002). First, we ordered the site-by-species incidence matrices, with sites a priori sorted by depth, with reciprocal averaging (RA) in correspondence analysis (CA). Then, objective criteria based on coherence, turnover, and boundary clumping were used to assess the correspondence of the empirical data set with each of the hypothetical idealizations of species distribution, that is checkerboard, nested, evenly spaced, Gleasonian, or Clementsian patterns (Leibold & Mikkelson, 2002;Valanko et al., 2015). The significance of the index value for coherence and turnover was tested using a fixedproportional null model, which maintains species richness of each site (i.e., row sums are fixed), but species ranges are filled based on their marginal probabilities (i.e., the "r1" null model, Gotelli, 2000;Dallas, 2014). We used 1,000 simulations to provide random matrices, with zero rows allowed in the null matrices because of sparse measured site-by-species matrices (Dallas, 2014). The observed index values were compared to the distribution of index values derived from randomization with the null model to assess statistical significance. We interpreted the results of the EMS analysis according to  and used the metacommunity function in the "metacom" package (version 1.5.0, 2018-01-20) for calculations and the "Image" function for graphs of sorted matrices (Dallas, 2014), in the R environment (R Core Team R version 3.4.3, 2017-11-30). The data used for the EMS analysis were binary, that is presence-absence of species. Pearson r correlation was used to test whether latent main environmental gradient (i.e., primary axis site scores from RA in correspondence analysis) was significantly correlated with the measured predictor variable (i.e., water depth).

| Test of potential bias due to different number of sites and different distributions of site depths
The 1880s dataset contained more sampling sites and a higher proportion of deep sites than the recent datasets (Supporting Information Table S1). To evaluate a potential bias of the EMS results and total richness, we therefore permuted station sets in the old dataset with the same number of sites (19) as in the recent datasets, and for each permutation constructed new site-by-species matrices subsequently subjected to EMS analysis in order to determine the depth structure pattern. We applied Pearson r correlation analysis between the site scores on the main RA axis and depth for each permutation. This allowed us to estimate the likelihood of the same pattern in the 1880s using 19 sites as in recent times. The permutations, 600 times, were performed both on the 1800 dataset with all species and the set with species occurring both in old and in recent times.

| Community changes
The total number of species from the three decades was similar: 239 in the 1880s, 236 in the 1990s, and 241 in the 2000s. Of these, only ca. half-117 species-were recorded both from the 1880s and from one or both of the two recent decades (Supporting Information   Table S2). Thus, species turnover between old and recent times was ca. 50%.
Using all species, there were highly significant coherent patterns in all three decades, that is fewer embedded absences than predicted by the null model (Table 1). Turnover, the number of replacements, was higher than predicted from the model but insignificant in the 1880s, while turnover was significantly higher than expected in the two recent decades ( Table 1). Clumping of ranges was significantly higher in both recent decades. Thus, we classified the structure in the 1880s as Quasi-nested, which means that most species had completely overlapping depth distributions with many intervals reaching over the whole studied depth range and that species depth ranges were close to randomly distributed along the main RA axis (Figure 4a). In contrast, we classified the states in the 1990s and the 2000s as Clementsian, which means that many species had narrower depth distribution intervals partly overlapping and with groups of species with similar ranges ordered along the main RA axis  Table 2). The only difference was that the 1880s dataset did not show significant coherence and was therefore indistinguishable from random.
While the site scores on the primary RA axis were virtually unrelated to water depth in the 1880s, the scores from the two recent patterns were highly significantly (p < 0.001) positively correlated with water depth ( Figure 5).

| Test of potential bias due to different number of sites and different distributions of site depths
To test the robustness of the difference in depth structure between old and recent times, we permuted station sets in the old dataset with the same number of sites (19)

| Change in depth distribution of individual species occurring both in old and in recent times
The maximal depth of occurrence, that is the lower depth limit of the depth range, decreased for 58 species in the 1990s and 50 species in the 2000s with on average ca. 20 m for most species up to, or above, 30 m depth (Table 4, Figure 6). This species category included species from all six main taxonomic groups, and more than half of the species belonged to Bivalvia, Gastropoda, and Decapoda.
Concurrently, 41 species in the 1990s and 36 in the 2000s, in the 1880s mainly found in the HFP zone, increased their maximal depth of occurrence with at least ca. 20 m, many to the greatest depth of the study (Table 4, Figure 7). The maximal depth limits of remaining species, 10 species in 1990s and 11 species in 2000s, were unchanged relative to the 1880s. While the decreasing group contained species from six higher groups, the increasing group contained species from four groups (no gastropods or decapod crustaceans at all), and 74% of the species belonged to Polychaeta and Peracarida crustaceans.

| Change with respect to traits with different sensitivities to bottom trawling
We prepared Supporting Information Table S3 with Table 4).

| Species richness
Richness is dependent on sampling effort, and because, for example, the number of sites differs greatly between the old and recent datasets, this precludes direct comparison between the datasets.
Therefore, using the rarefaction in the above-described permutation test, we estimated richness in 19-site samples in the old datasets for comparison with the recent datasets (Table 5). Richness in the 1880s estimated from 19-site samples was considerably lower than measured richness in the recent datasets, which were well above the upper limit of the estimated mean + 1SD (Table 5). Much of this difference, however, was due to a higher fraction of polychaete species in recent datasets, 44%-48% compared to 18% in the 1880s.
Assuming that the same fractions of polychaetes and nonpolychaetes in the estimated mean richness as in the measured richness from 1880s indicate similar level of richness of the nonpolychaete fauna in old and recent times, richness of polychaetes increased threefold in recent times (Table 5).

| D ISCUSS I ON
The different sampling gears in old and recent times makes quantitative comparisons between periods problematic. The dredges and trawls used in the 1880s likely sampled large, epifaunal organisms with low density much more efficiently than the grabs used in recent times, while the grabs used in recent times probably sampled abun- showing increasing maximal depth between old and recent times. More information in legend of Figure 6 Part of the high richness of polychaetes in recent times may also be due to taxonomic evolution. Species from several genera like Nepthys, Pholoe, Chaetopterus has been split into more species rather than merged into fewer species, and several new polychaete species have been described later than the surveys in 1880s (e.g., Hartmann-Schröder, 1996). To overcome some of these methodological biases, we, in addition to using binary data only, also compared distributions between time periods of those species actually recorded both in old and in recent times. For the depth distribution structure of all species in each dataset however, we believe that there is little reason to assume a bias as long as sampling methodology was consistent within each period.
The depth distribution structure has changed from a Quasinested/random structure in the 1880s with great overlaps be-  in the adjacent Baltic Sea (Sanden & Håkansson, 1996;ca. 5 m in 100 years). The depth penetration of the autotrophic Zostera marina in the Kattegat has decreased from at least 20 m in the 1880s (Petersen, 1893) to 10 m or less in recent times (Duarte, 1991;Nielsen, Sand-Jensen, Borum, & Geertz-Hansen, 2002).
Although the group of species with shallower distribution in recent times contained most large bivalves (Supporting Information   (Conley et al., 2007;Rosenberg, Loo, & Möller, 1992). The 1990s dataset was sampled 10 years after the main hypoxic event in 1988 (Rosenberg et al., 1992), and the 2000s dataset 6 years after the 2002 event (Conley et al., 2007). Recovery time after anoxia-induced extermination of the fauna in some Swedish coastal areas may be some 5-6 years for most species (e.g., Josefson, Blomqvist, Hansen, Rosenberg, & Rygg, 2009 (Figure 3) with peak intensity in fair agreement with Pommer et al. (2016). It may therefore be more than a coincidence that the maximal depth occurrence of many species has moved upward to, or above 30 m, away from the depth zone intensively trawled today.
We find that the composition in terms of traits differs markedly between most species with shallower occurrence than 30 m in recent times and species that have remained within the depth zone deeper than 30 m (Figure 8). A majority of the species in the first mentioned group (green in Figure 8) are emergent or semiemergent and contain most species with large maximal body size in this trait category, traits that make them sensitive to trawling (Kaiser et al., 2000; Supporting Information Table S3 our study).
In contrast, the infauna trait group is highly dominated by species remaining in the deep zone fished today (brown in Figure 8), particularly small-sized infauna. The differences in trait composition between the species with shallower occurrence and species remaining in the fished zone reflect to some extent differences in taxonomic composition between decreasing and increasing species with respect to maximal depth. Large-bodied emergent/ semiemergent species contain most large bivalves, gastropods, and decapods, while small infauna mostly are polychaetes. Of the species category with deeper distribution in recent times, but still occurring in the HFP zone, a majority belonged to Peracarida crustaceans and polychaetes. Most species of these groups are small and short-lived, many polychaetes with a longevity of ca. 2 years and several Peracarida crustaceans like the Cumacea species even annual. Polychaetes are truly infaunal although several feed at the sediment-water interface. We expect species with these traits to survive in trawled areas (e.g., Kaiser et al., 2000). For instance, the fishing gears are less likely to injure small-bodied species because of displacement by the pressure wave in front of the gear (Gilkinson, Paulin, Hurley, & Schwinghamer, 1998), and smallbodied taxa are typically able to recover fast and hence withstand higher levels of chronic trawling.
All Peracarida crustaceans "moving" from the fished zone to the unfished zone in one or both of the recent decades were semiemergent and included several species making tubes or pillars protruding above the sediment surface, that is from the genera Dyopedos, Ampelisca, and Haploops. The last mentioned formed extensive tube mats in deeper areas in the pretrawling period (Petersen, 1913) and  Table   S2). The disappearance of Haploops tube mats in the Kattegat is probably one obvious effect of demersal trawling, because mats, although with local extent, still existed in the 2000s in the Sound (Josefson & Pedersen, 2001) with a similar physical environment but with a trawling ban.
So, the decreases in maximal depth away from the fished zone, and the different species and trait composition of decreasing and species remaining in the fished zone are as argued above, are most easily explained by benthic trawling impacts in recent times.
Concurrently, richness in the investigation area corrected for sampling effort did not decrease, and therefore, it seems that the major effect of trawling is changes in spatial distribution of species within the area rather than net species loss.
To conclude, the depth distribution structure of benthic invertebrate fauna in the eastern Kattegat has changed substantially since the 1880s. Many species, including the species likely to suffer most from intensive trawling, have narrower depth distribution intervals in recent times with changed maximal depth limits above and away from the fished zone. Our study using a temporal preimpact reference underpins earlier suggestions that chronic fishing disturbance has changed benthic community structure over wide areas of the shelf seas (e.g., Kaiser et al., 2000). We suggest that the Petersen dataset from the 1880s, despite some taxonomic and methodological caveats, is highly suitable for setting reference conditions for management of benthic species distributions in this sea area.

ACK N OWLED G M ENTS
We thank current and past employees at the Swedish Environment Protection Agency and the Swedish Counties for sampling and determination of samples from recent times, which has made the present comparison possible. Thanks go to M. Jahnke and two anonymous reviewers for valuable comments on the manuscript, which greatly improved the manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S CO NTR I B UTI O N
AJ and LL gave the idea and conceived the study; AJ wrote the paper; LL, MB, and JR collated the data; AJ and MB analyzed the data. All performed critical reading and contribution to the revisions.

DATA ACCE SS I B I LIT Y
All primary station and fauna data are contained within the paper and its appendix files.