Marine long-term biodiversity assessment suggests loss of rare species in the Skagerrak and Kattegat region

Studies of cumulative and long-term effects of human activities in the ocean are essential for developing realistic conservation targets. Here, we report the results of a recent national marine biodiversity inventory along the Swedish West coast between 2004 and 2009. The expedition revisited many historical localities that have been sampled with the same methods in the early twentieth century. We generated comparable datasets from our own investigation and the historical data to compare species richness, abundance, and geographic distribution of diversity. Our analysis indicates that the benthic ecosystems in the region have lost a large part of its original species richness over the last seven decades. We find evidence that especially rare species have disappeared. This process has caused a more homogenized community structure in the region and diminished historical biodiversity hotspots. We argue that the contemporary lack of rare species in the benthic ecosystems of the Kattegat and Skagerrak offers less opportunity to respond to environmental perturbations in the future and suggest improving the poor representation of rare species in the region. The study shows the value of biodiversity inventories as well as natural history collections in investigations of accumulated effects of anthropogenic activities and for re-establishing species-rich, productive, and resilient ecosystems.


Introduction
Marine habitats experience rapid declines in biodiversity worldwide (Jackson et al. 2001;Halpern et al. 2008), a process that creates urgent demand for a better understanding of the long-term effects of such severe alterations in ecosystem diversity (Rockström et al. 2009;Lotze 2010). Over the past decades, the impacts of major anthropogenic pressures on coastal and benthic marine biodiversity have been studied intensely. Here, especially field assessments investigated the negative effects arising from bottom trawling (Jennings and Kaiser 1998;Kaiser et al. 2006;Tillin et al. 2006;Worm et al. 2006;Olsgard et al. 2008), coastal nutrient loading (Rosenberg and Nilsson 2005;Quijon et al. 2008), and climate change (Norderhaug et al. 2015). However, since many of these drivers act on the ecosystems simultaneously and over long periods of time Halpern et al. 2008), it is difficult to infer cumulative impacts from such assessments (Moksnes et al. 2008;Robinson and Frid 2008). Furthermore, in many experimental field studies, the impacts are already a Communicated by D. M. Paterson Electronic supplementary material The online version of this article (doi:10.1007/s12526-017-0749-5) contains supplementary material, which is available to authorized users. part of the control (Pauly 1995), and hence it is not possible to rely on contemporary and experimental investigations alone when examining long-term changes in marine ecosystems.
Historical studies can offer valuable insight in ecosystemwide responses to the overall sum of human pressures that act in concert and over long periods of time. For such investigations the benthos around the region of the North Sea is especially well suited. Compared to most coastal regions of the world, this area has a long history of biological recording with quantitative surveys dating back to the late nineteenth and early twentieth century (Petersen 1918;Robinson and Frid 2008;Narayanaswamy et al. 2010). Based on this information, a large number of historical comparisons have already been carried out. For example, long-term investigations of fish, plankton, and benthos in the Western English channel found indications for regime shifts during the last century caused by fishing pressures (Southward et al. 2005). Other studies investigated the Southern and central parts of the North Sea and found similar evidence for long-term changes in benthic community structure attributed to fishing (Pennington et al. 1998;Rumohr and Kujawski 2000;Bradshaw et al. 2002;Robinson and Frid 2008) and nutrient loading (Schroeder 2005;Schumacher et al. 2014). Likewise, a series of historical studies from the eastern parts of the North Sea documented remarkable long-term changes in benthic communities attributed to trawling pressure and eutrophication (Rosenberg and Möller 1979;Pearson et al. 1985;Rosenberg et al. 1987;Göransson 2002). Such assessments allow valuable insight into long-term transformation of ecosystems and may help to identify important biodiversity trends in a region, a feature highly relevant for future conservation policies (Pereira et al. 2013).
During the early twentieth century, an expedition led by L.A. Jägerskiöld inventoried the benthic diversity of the Kattegat and Skagerrak region (Jägerskiöld 1971). We were able to revisit many of these historical locations during a recent national marine biodiversity inventory program, and we used this opportunity to test the validity of historical data for understanding the long-term trends in biodiversity in this region. Our expedition re-sampled a large number of historical locations with similar equipment and methods, and we subsequently generated comparable datasets from both expeditions to analyse ecosystem-wide changes in species diversity over a period of more than 70 years.

Data collection
The historical inventory was carried out in the Swedish, Danish, and Norwegian Economic Zone by L. A. Jägerskiöld in 1921-1938(Jägerskiöld 1971. Overall, 440 benthic localities were visited in the Kattegat and Skagerrak, from shallow to deep water, and typically between spring and autumn ( Fig.  1a), generating a dataset with 33,661 species observations. Usually, several samples were taken per locality by combining various dredges and trawls (Tables 1 and 2). Organisms that were obviously picked up in the water column were discarded. All living marine invertebrate species larger than 1 mm were collected, identified, preserved, and stored (fixed in formaldehyde, ethanol, or formol, and later transferred to ethanol). Subsequently, all specimens were vouchered, catalogued, reexamined by experts if necessary, and finally stored at the Gothenburg Natural History Museum, Sweden.
Our recent inventory revisited the historical locations between 2004 and 2009 (Fig. 1a). Here, altogether 504 localities were sampled during spring and autumn, generating a dataset with 17,249 species observation records. The equipment used was of the same type and with similar dimensions and mesh size as applied in the historical inventory (Tables 1 and 2). The same criteria for collection and identification were applied as in the historical inventory, while all material was stored in ethanol at the Gothenburg Natural History Museum, Sweden.
We prepared and submitted all original data from our own inventories as well as from the museum collections ( Fig. 1a) including overall 50,910 species records to the Swedish Environmental and Climate Data Repository (www.ecds.se) under the identifier 01d148f8-f87c-47e3-adfc-c10619f6e9a1 (metadata) and http://webdav.swestore.se/snic/ecds/prod/ BenthicInventories/ (data file). Data cleaning, refinement, and taxonomic name resolution We employed a semi-automated workflow developed by Mathew et al. (2014) to generate comparable datasets from the metadata of both inventory programs. The workflow provides a user interface for preparation of taxonomically accurate species lists and observational records and can be executed online (https://portal.biovel.eu/workflows/641). All locations were assigned to the following habitat categories: soft bottom (sand, mud), hard bottom (rocks, boulders, stones), or shell gravel. We then structured species occurrences and sample locations according to habitat, geographical reference, sampling gear, and depth profile. Revisited and comparable samples were defined as locations with similar geographical reference up to the first decimal of both latitude and longitude. Locations also had to have overlapping or adjacent depth profiles, while all locations above the halocline (app. 15 m) were excluded because the research vessels could not adequately sample such shallow habitats. Based on these criteria we selected a group of 54 revisited and comparable localities in the Swedish Exclusive Economic Zone (Fig. 1b, Table 2).
Sample effort was calculated as overall haul volume from the dimensions of the individual sample equipment and the haul length at each location (Tables 1 and 2). For the historical samples, this information was derived from the original logbooks of the Jägerskiöld inventory available at the Gothenburg Natural History Museum, Sweden.
Species observations from selected localities were cleaned and refined in the following order. We only included taxa unambiguously identified by taxonomic experts in both inventories, but excluded endoparasites. Spelling errors and variations in species names were identified and corrected with the taxonomic 'cluster' function of the workflow (Mathew et al. 2014). A few ambiguous entries such as missing species epithet (only listed as sp.), records with genus names only, or entries with 'cf' references were either resolved or excluded. All species names were transformed to the accepted name provided by the web services of Catalogue of Life (Roskov et al. 2014) and World Register of Marine Species (WoRMS) to eliminate all synonymies in the dataset (Mathew et al. 2014). Multiple records of the same species in a location were removed, leaving only presence information for further analysis in the data set. The resulting dataset (shown in Fig. 1b, available as supplementary online material File S1) was then used for the statistical analysis.

Statistical analysis
The variation between the samples was explored to identify and correct for any potential sampling bias in the historical and recent inventories. Variables included in the analysis were species richness, sampling effort, habitat, geographic location, depth, and season. First, we plotted the frequency distribution of mismatches in sampling effort, average sampling depth, sampling date, and substrate type for all revisited localities (Fig. 2). Second, we performed multidimensional scaling of all samples to visualize the similarity in species composition among and between localities of the two inventories (Fig. 3). Finally, we log transformed species richness and sampling effort (as advised by a preliminary Box-Cox analysis) and evaluated the variables using a generalized additive model (GAM) with bi-dimensional smothers, taking the dishomogeneity in the sampling effort into account ). Altogether, 703 models (160 of which included a biodiversity hotspot variable defined after exploratory analysis) were tested to compare the effect of the different variables on the variance in the data set ( Fig. S1, Table S1-S3).
Species richness (number of species) was calculated as a function of sample effort (haul volume). In addition, we applied a classical rarefaction analysis, re-sampling the observations and the localities to assess significance and the species richness abundance-based coverage estimator (ACE) for each inventory. We used a non-parametric species richness  Table 2 Overview of selected localities and sample properties in the refined dataset. Column legend: 1, locality station code; 2, geographic coordinates; 3, minimum-maximum depth (in m); 4, date; 5, habitat type (S, soft bottom; H, hard bottom; G, shell gravel); 6, sample equipment (with haul length in m), see Table 1 for specifications of individual sample volume; 7, sample effort as overall haul volume per location (in m 3 ); 8, species richness (no. of species); 9, relative species richness (no. of species/100 m 3 haul volume); 10, distance between localities (in m) Historical inventory  estimator (ACE-1) as advised in Gotelli and Colwell (2011) to obtain a quantitative estimate of species richness for repeated incidence data. Analyses were implemented with R scripts (Wang 2011;R Core Team 2013). Species abundance was calculated as the relative frequency of occurrence for individual species in the selected 54 locations. We defined abundance thresholds at 10% of the total number of locations for rare species and at 50% for common species. This allowed us to classify all species into four abundances classes: absent (0%), rare (0% -<10%), intermediate (10-50%), and common (>50%).
Geographical structure was assessed by comparing the marginal values of the Akaike information criterion weight (wAIC) across 543 generalized additive models resulting from the different parameterisations of the geographical information, including spatial smothers (Table S1). Subsequently, the geographical structure was modelled with groupings of historical hotspots for a total of 160 models (Table S2). Historical hotspots were defined as the locations with the highest species richness in the historical data set, and included 16 localities with an average richness of 103 species/100 m 3 haul volume. The analysis was done with R scripts (Wood 2006;Rhodes et al. 2009; R Core Team 2013).

Exploration of variance
Variation in sample effort was considerable (Figs. 2a, 4b, Table 2), both within and between inventories (historical/recent effort minimum = 23/16 m 3 , maximum = 322/322 m 3 , average = 104/78 m 3 , standard deviation = 80/60 m 3 ). Most of the revisited localities had a mismatch in sampling effort between 20 and 100 m 3 (Fig. 2a). Depth profiles only had substantial variation within each inventory, but not between inventories (historical/recent depth minimum = 11/19 m, maximum = 200/220 m, average = 51/47 m, standard deviation = 36/35 m). Most of the revisited localities showed mismatches in depth of ±15 m (Fig.  2b). The seasonal variation within and between inventories was confined to the summer months. Historical locations were sampled mostly between June and August, while most of the recent localities were sampled either in May or in August, resulting in a seasonal mismatch between the revisited localities of up to 3 months (Fig. 2c). The variation in the substrate between historical and recent localities was very small. Most localities consisted of either soft bottom or mixed soft and hard bottoms at both times of sampling (Figs. 2d and 3). Overall, there was no bias in the data sets that indicates systematic over-or under-sampling of a certain habitat, season, depth, effort, or sample gear in one of the two inventories (Figs. 2 and 3) and localities are comparable with regard to these variables. The variation in sampling effort, however, needs to be included when comparing species richness. The GAM approach attributed the variance to the habitat, biodiversity hotspots, sampling effort, and depth. Apart from depth, all variables had a different influence between the historical and recent inventory (Fig. S1, Table S3). In the historical inventory, 16.2% of the variance in species richness was assigned to sampling effort (significant at P = 0.00277), while in the recent inventory, there was no relationship between species richness and sampling effort (P = 0.95). This indicates that the sampling in the historical inventory was still unsaturated (Fig. 4b, Table S3).

Species richness
The refined dataset included, overall, 648 species (4412 species records) across 54 revisited locations in the investigated region ( Table 2). The historical partition included 607 species (3282 species records), while the recent partition included 254 species (1130 species records). Overall, 32.8% of the species were recovered across both inventories, while 60.8% occurred exclusively in the historical dataset, and 6.3% occurred only in the recent dataset.
The rarefaction curve indicated a halfway reduction of recent compared to historical species richness (Fig. 4a), which was confirmed by the ACE-1 estimator (Chao and Lee 1992) showing that the overall estimated species richness in the recent data set decreased to 47.2% of the historical values (Fig. 5). On average, historical localities recovered 88.2 species/100 m 3 haul volume, while recent localities recovered only 38.6 species/100 m 3 haul volume (Table 2). Also, the relations between species richness and sample effort were different between the historical and recent data set. Increasing effort resulted in higher species richness in historical samples, but this was not the case in recent samples (Fig. 4b).

Abundance trends
The comparison of abundances between the historical and recent data sets showed a predominating negative trend (Fig. 6, Table 3). Overall, 74.7% of the investigated species  Among the species with strongest positive abundance t r e nd s , w e fo un d c or a l s (C a r y o p h y l l i a s m i t h i i , Kophobelemnon stelliferum, Alcyonium digitatum), bivalves (Nucula nitidosa N. nucleus, Mysia undata, Abra alba, Pecten maximus, Thracia convexa, Pododesmus patelliformis, Modiolarca subpicta), and one polychaete (Nephtys kersivalensis). These species are typically small-to mediumsized suspension or deposit feeders, living on top or within the sediment (epifauna, infauna). Among the species with the strongest negative abundance trends, we found especially polychaetes (Pectinaria auricoma, Goniada maculata, Aphrodita aculeata, Pista cristata, Owenia fusiformis), but also one echinoderm (Psammechinus miliaris), a crustacean (Verruca stroemia), and a mollusc (Buccinum undatum). These species are typically small-to medium-sized predators, scavengers or suspension feeders, living on top or within the sediment (epifauna, infauna).

Geographic distribution of diversity
The analysis of the historical samples across 160 GAM models with a biodiversity hotspot predictor rendered three areas with high biodiversity (n = 16) in the historical data set (Fig. 7). The northernmost area (n = 7) lies in the centre of a national marine sanctuary (Gonzalez-Mirelis et al. 2014), while the others are located on the shallow water banks (n = 3) and in the coastal zones (n = 6) between Denmark and Sweden. By grouping all 16 species rich localities together, we obtained a 95% confidence set with 6 models, while the best model included 54% of the weight of the selected models (Table S2). The geographical structure is well described by the 2172 Mar Biodiv (2018)

Discussion
Sampling bias and working with inconsistent data sets Our study investigated the changes in richness, abundance, and geographic distribution of benthic species in the Kattegat and Skagerrak region based on the recordings made by two large biodiversity inventories. Although both inventories had a very similar design, the data generated by these investigations retained considerable heterogeneity. By using a carefully filtered subset of 54 revisited localities sampled with similar methods during the summer season, we could remove a large part of the heterogeneity. However, some variance in depth, season substrate, and constitution of the sample gear remained and may influence the observed patterns. However, we found no evidence of systematic bias in any of these variables across both inventories. Additional factors that may have influenced the historical changes observed in this study may be the coarse partitioning of the localities into three habitat types (hard bottom, soft bottom, and shell gravel), or the estimation of sample effort from the equipment's catchment area and the overall haul length. Also, the overall sampling period, which was longer for the historical data set  than for the recent one (2004-09) may have captured more of the temporal species turnover and hence contributed to the higher levels of alpha diversity observed in the historical data. In conclusion, our results cannot entirely be assigned to the long-term changes in the region, and may to some extent be caused by the deviations between the sampling routines that we could not control. In this context, our objective was not to remove all sampling bias from the data, but attempt to find a trade-off between the degree of harmonization of the data sets and the conclusions that can be drawn from the analysis.

Historical trends
Our results indicate a reduction of species richness in the investigated region, which can be explained to some extent by the extirpation of rare species. The comparison of species abundances in the historical and recent inventory suggests that more than 50% of the recorded species change status from rare to absent. This trend seems to have contributed to a more homogenized community structure across the investigated region and may have diminished former biodiversity hotspots. Although the longer sampling period of the historical inventory may have increased the overall capture of rare species, it does not influence the species richness at the level of individual localities. Our results show that historical localities not only have a higher yield of species for a given sampling effort, but also a tendency to obtain more species with higher sampling effort. This is not the case in recent localities, where higher sampling effort does not yield more species. One of the explanations for these differences may be the absence of rare species in the recent localities. A comparison with similar historical accounts shows that substantial regime shifts over longer time periods are known for the region (Rosenberg and Möller 1979;Robinson and Frid 2008). Our results show a 32% overlap in species assemblages between the historical and recent inventory, a figure that is very similar to estimates obtained by Pearson et al. (1985) and Rosenberg et al. (1987) over a similar period. However, in contrast to previous long-term studies, which report the species turnover as a balance between species recruitment and extirpation, our study also indicates a severe loss of alpha diversity in addition to the regime shift.

Potential causes and consequences of species reduction
A deeper analysis of ecological responses due to changes in the species composition is limited by the lack of consistent trait information for the majority of the species included in this study. However, some responses to prevalent pressures in the region may be discussed using representative species. Anthropogenic drivers that influence benthic diversity in the North Sea are typically associated with pollution, overfishing, habitat destruction, invasive species, and climate change Doney and Schimel 2007;Halpern et al. 2008). The reduction of alpha diversity observed in our study is likely to be related to a combination of these factors, but we have little indication of any specific causes for the observed loss in species richness. However, depletion of ecological niches through continuous physical disturbance (e.g. seabed trawling) is a well-documented process in the region, which is known to remove rare and specialist species from the ecosystem (Jennings and Kaiser 1998;Kaiser et al. 2006;Clavel et al. 2011). In addition, we observed that fragile infauna such as burrowing and tube-dwelling polychaetes were among the species with the strongest signals of decline, which may also be a response to continuous trawling activities. In contrast, more robust in-and epifauna like molluscs and corals seemed to have increased over time, and this pattern may be indicative for a better resistance of such species to physical disturbance. However, the observed changes in community composition are likely to be caused by a combination of several factors. Many of the species that have become more dominant are suspension and deposit feeders, organisms that benefit from the elevated nutrient levels in the region (Graneli and Sundback 1985;Posey et al. 1999;Karlson et al. 2002), while other important drivers of the observed changes in community composition may also include constantly rising sea temperatures and associated changes in water quality that can affect individual species differently (Hiddink et al. 2015). In conclusion, our study only includes two distant sampling periods, while specific environmental data associated with potential drivers of change were not available for analysis. Future inventories and monitoring programs should therefore emphasize the collection or linkage to environmental information to allow a deeper understanding of the impact of specific drivers on biodiversity decline.
The reduction of rare species indicated in our study may have implications on the resilience of the benthic ecosystems in the region. Specialist decline is known to cause functional homogenisation, which effects ecosystem functioning and productivity in the long-term and may ultimately lead to deterioration of important ecosystem services (Clavel et al. 2011;Cardinale et al. 2012). An increasing amount of studies shows that undisturbed marine ecosystems possess a broader functional reservoir allowing them to react better to environmental perturbations compared to exploited systems (Stachowicz et al. 2007;Rasher et al. 2013). Therefore, the constitution of rare species in the ecosystem should be better monitored in marine conservation programs. Appropriate indicators for biodiversity already exist as a part of national and international legislation (Borja 2006) and should be used to actively improve the representation of rare species in benthic assemblages. This intention could be realized by introducing new assessment methods into biodiversity monitoring programs, which enable a better accounting for rare species (Bourlat et al. 2013).