Modelling and predicting habitats for the neobiotic American razor clam Ensis leei in the Wadden Sea

Journal Pre-proof Modelling and predicting habitats for the neobiotic American razor clam Ensis leei in the Wadden Sea Philipp Schwemmer, Sven Adler, Leonie Enners, Henning Volmer, Johanna Kottsieper, Klaus Ricklefs, Maria Stage, Klaus Schwarzer, Kerstin Wittbrodt, HansChristian Reimers, Kirsten Binder, Ragnhild Asmus, Harald Asmus, Sabine Horn, Ulrike Schückel, Jörn Kohlus, Kai Eskildsen, Knut Klingbeil, Ulf Gräwe, Stefan Garthe PII:

S0272-7714(19)30191-X

DOI:

https://doi.org/10.1016/j.ecss.2019.106440

Reference:

YECSS 106440

To appear in:

Estuarine, Coastal and Shelf Science

Received Date: 27 February 2019 Revised Date:

2 September 2019

Accepted Date: 17 October 2019

Please cite this article as: Schwemmer, P., Adler, S., Enners, L., Volmer, H., Kottsieper, J., Ricklefs, K., Stage, M., Schwarzer, K., Wittbrodt, K., Reimers, H.-C., Binder, K., Asmus, R., Asmus, H., Horn, S., Schückel, U., Kohlus, Jö., Eskildsen, K., Klingbeil, K., Gräwe, U., Garthe, S., Modelling and predicting habitats for the neobiotic American razor clam Ensis leei in the Wadden Sea, Estuarine, Coastal and Shelf Science (2019), doi: https://doi.org/10.1016/j.ecss.2019.106440. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

1 1

Modelling and predicting habitats for the neobiotic American

2

razor clam Ensis leei in the Wadden Sea

3 4

Philipp Schwemmer1, Sven Adler2, Leonie Enners1, Henning Volmer1, Johanna

5

Kottsieper1, Klaus Ricklefs1, Maria Stage1, Klaus Schwarzer3, Kerstin Wittbrodt3,

6

Hans-Christian Reimers4, Kirsten Binder4, Ragnhild Asmus5, Harald Asmus5,

7

Sabine Horn5, Ulrike Schückel6, Jörn Kohlus6, Kai Eskildsen6, Knut Klingbeil7, Ulf

8

Gräwe7, Stefan Garthe1

9 10

1

11

University of Kiel, Hafentörn 1, 25761 Büsum, Germany

12

2

13

Management, Division of Landscape Analysis, 901 83 Umeå, Sweden

14

3

15

of Kiel, Otto Hahn Platz 1, 24118 Kiel, Germany

16

4

17

Hamburger Chaussee 25, Flintbek

18

5

19

Wadden Sea Station Sylt, Hafenstraße 43, 25995, List, Germany

20

6

21

Conservation – National Park Authority, Schlossgarten 1, 25832 Tönning

22

7

23

Baltic Sea Research Warnemünde (IOW), Seestraße 15, 18119 Rostock, Germany

24 25

Research and Technology Centre (Forschungs- und Technologiezentrum),

Swedish University of Agricultural Sciences, Department of Forest Resource

Institute for Geoscience, Sedimentology – Coastal and Shelf Research, University

State Agency for Agriculture, Environment and Rural Areas of Schleswig-Holstein,

Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung,

Schleswig-Holstein Agency for Coastal Defence, National Park and Marine

Department for Physical Oceanography and Instrumentation, Leibniz Institute for

2 26

Corresponding author:

27

Philipp Schwemmer, Research and Technology Centre (Forschungs- und

28

Technologiezentrum), University of Kiel, Hafentörn 1, 25761 Büsum, Germany

29

Fax: +49 4834 604 299; Email: [email protected]

30 31

Abstract

32

Neobiotic species can have profound impacts on food webs and entire ecosystems.

33

The American razor clam Ensis leei was introduced into the Wadden Sea by

34

vessels in the late 1970s and has since spread widely. It has been suggested that

35

Ensis does not interact strongly with other benthic species. The abundance and

36

biomass of E. leei were recorded in 2,393 samples in the north-eastern Wadden

37

Sea and 800 samples in the south-eastern Wadden Sea over a total period of 9

38

years. Using an interdisciplinary approach, we developed a habitat prediction

39

model using sedimentological and hydrodynamic predictors to help understand the

40

shape of the ecological niche occupied by Ensis in the Wadden Sea. Our model

41

showed that Ensis preferred areas with moderately high bed shear stress and

42

prolonged or constant water coverage. Ensis preferred coarse sediments in the

43

northern sub-area but coarse and muddy sediments in the southern sub-area and

44

was negatively affected by the sand mason worm Lanice conchilega in the northern

45

sub-area. Predictions of the spatial distribution of Ensis using the northern and

46

southern datasets revealed major differences in predicted hot-spots throughout the

47

entire study site. This study thus highlights the need to collect a sufficiently large

48

dataset from different sub-areas of the Wadden Sea to allow valid conclusions to

49

be drawn regarding the spatial distribution of Ensis. The negative effects of L.

50

conchilega on Ensis abundance and biomass as well as the occurrence of Ensis in

3 51

muddy sediments in the south suggest that the ecological niche of this neobiotic

52

species is likely to overlap partly with the native fauna of the Wadden Sea.

53 54

Key Words: habitat model; invasive species; ecological niche; hydrodynamics;

55

sediment; Lanice conchilega

56

4 57

1. Introduction

58

An increasing number of neobiotic species has been recorded in different marine

59

systems throughout central Europe (e.g., Sahlin et al., 2011; Galil et al., 2014). This

60

is also true for the World Heritage Site of the Wadden Sea (e.g. Buschbaum et al.,

61

2012), which is a relatively young ecosystem with a consequently low species

62

diversity and high environmental variability facilitating a broad range of ecological

63

niches (Reise, 1985). However, invasive species have the potential to alter food

64

web dynamics, biodiversity and features of the seafloor in the Wadden Sea (Baird

65

et al., 2012; Lackschewitz et al., 2015). Although invasive species in the Wadden

66

Sea have not yet been described to cause profound damage to the ecosystem,

67

future changes under an altered climatic or hydrodynamic regime cannot be ruled

68

out (Reise et al., 2006). This emphasises the necessity to explore the ecological

69

niche of invasive species in the new system they occupy.

70

The American razor clam Ensis leei (formerly Ensis directus, hereafter referred to

71

as Ensis) is rated as a widespread species in Western European regions

72

(Severijns, 2002; Galil et al., 2014). It was introduced into the Wadden Sea by

73

vessels and its widespread distribution was first documented in the early 1980s

74

(Armonies & Reise, 1999; Dekker and Beukema, 2012). Further expansion of its

75

distribution throughout European seas is highly likely and favoured by climate

76

change (Raybaud et al., 2015). Due to its high biomass in the Wadden Sea and its

77

high potential to expand, Ensis has been rated as an invasive species with

78

potentially strong impacts (Gollasch and Nehring, 2006). The aim of this study is

79

therefore to develop a model on the habitat potential of this neobiotic species which

80

can be used to inform the descriptors of invasive species (D2), seafloor integrity

81

(D6) and biodiversity (D1), as well as food webs (D4), within the Marine Strategy

5 82

Framework Directive (MSFD). This should be achieved by developing a prediction

83

model for Ensis, valid throughout the eastern Wadden Sea and adjacent offshore

84

areas with the potential also to transfer it to the entire Wadden Sea. This predictive

85

model could help to identify the ecological niche occupied by Ensis by quantifying

86

thresholds of environmental predictors that determine its occurrence. This will

87

contribute to understanding the role of Ensis in the whole food web (D1 and D4).

88

Furthermore, we discuss how strongly this niche might overlap with the native

89

fauna of the Wadden Sea (D2 and D4) and how this species is related to certain

90

properties of the seafloor (D6).

91

Ensis has its native distribution in the north-west Atlantic from the Gulf of St.

92

Lawrence along the whole US coast to the Gulf of Mexico where it occupies mainly

93

shallow, sandy subtidal areas down to a water depth of 37 m (Theroux and Wigley,

94

1983; Leavitt, 2010). It has never been described to occur in muddy sediments in

95

its native range (Theroux and Wigley, 1983). Ensis prefers moderate water flow

96

velocities and reaches densities between 200 and 2,000 individuals / m² in its

97

native habitat (Leavitt, 2010). Shorebirds as well as snails and worms have been

98

reported to predate on Ensis in the west Atlantic (Leavitt, 2010). Ensis has

99

successfully invaded the Wadden Sea and adjacent offshore regions of the North

100

Sea (Gollasch et al., 2015), as indicated by its remarkably high biomass and

101

abundance (Tulp et al., 2010; Dannheim and Rumohr, 2012; Dekker and Beukema,

102

2012; Breine et al., 2018). Specimen in the Wadden Sea grow significantly shorter

103

than in its native range (von Cosel, 2009). It has been suggested that this species

104

does not interact strongly with other benthic species inhabiting the sublittoral area

105

of the eastern North Sea (Swennen et al., 1985; Dannheim and Rumohr, 2012;

106

Dekker and Beukema, 2012). However, recent studies found that razor clams were

6 107

an increasingly important food source for Wadden Sea birds such as the Eurasian

108

oystercatcher (Haematopus ostralegus; Swennen et al., 1985; Schwemmer et al.,

109

2016), herring gull (Larus argentatus, Cadée, 2000; Enners et al., 2018), and

110

seaducks (Tulp et al., 2010; Schwemmer et al., this issue), as well as for certain

111

fish species (Tulp et al., 2010).

112

Previous studies suggested that Ensis is one of the few benthic species that prefers

113

to inhabit coarse, mobile sediments, thus reducing competition with other bivalve

114

species (e.g., Armonies and Reise, 1999; Dannheim and Rumohr, 2012; Dekker

115

and Beukema, 2012; Gollasch et al., 2015). In addition to these empirical field

116

studies, experiments have also shown that razor clams seem to avoid areas with

117

low sediment redistribution (van der Heide et al., 2014). However, precise

118

information on the habitat requirements and habitat potential (i.e. areas most

119

suitable for settling) of this invasive species in the Wadden Sea is still scarce. The

120

current study therefore aimed to develop a habitat model to predict the abundance

121

and biomass of Ensis in the eastern part of the Wadden Sea and the adjacent

122

offshore zone using sedimentological and hydrodynamic variables. Similar

123

predictors have previously been used to forecast the species composition and

124

distribution of the macrobenthos in the Wadden Sea (Puls et al., 2012; Folmer et

125

al., 2016; Singer et al., 2016).

126

Our study focused on one site in the north-eastern part of the German Wadden

127

Sea, which is characterised by large tidal inlets, barrier islands, and the Halligen

128

(i.e. inhabited marsh islands that may be inundated during storm surges), and

129

another site in the south-eastern Wadden Sea, which is mainly shaped by the

130

influence of the Elbe Estuary causing a different morphology compared with the

131

north-eastern site (Fig. 1). We deliberately have selected two contrasting areas

7 132

which represent common features in the Wadden Sea (i.e. Barrier islands with

133

large tidal inlets as can be also found in East and West Frisia and river estuaries

134

such as the Weser and Ems estuary). Therefore, results of our study should be

135

transferrable to many parts of the entire Wadden Sea.

136

We proposed four hypotheses. 1) Given that Ensis has been reported to prefer to

137

live in coarse sands (e.g., Dannheim and Rumohr, 2012; van der Heide et al.,

138

2014),

139

abundance/biomass of Ensis and the proportion of fine sand and mud. 2) Ensis has

140

been reported to occur mainly in shallow subtidal and low-intertidal areas

141

(Armonies and Reise, 1999; Dannheim and Rumohr, 2012; Dekker and Beukema,

142

2012; Breine et al., 2018); we therefore expected to find a negative correlation

143

between Ensis and water depth as well as preference for subtidal areas and

144

intertidal areas showing long periods of inundation. 3) Ensis was shown to tolerate

145

mobile sediments (see above), and we therefore expected habitat choice to be

146

affected by bed shear stress (expressed as absolute force and duration of critical

147

bed shear stress; Kösters and Winter, 2014). Specifically, we expected to find an

148

optimum curve in Ensis abundance/biomass with respect to bottom shear stress,

149

given that areas with very low shear stress might not be mobile enough (compare

150

van der Heide et al., 2014), while very high shear stress might have critical effects

151

on sediment. 4) Despite morphological differences throughout the study area, we

152

assumed that correlations with environmental predictors would remain stable

153

throughout the area, enabling valid predictions to be made regarding the potential

154

distribution of Ensis. We tested this by setting up a prediction model using data

155

from the north-eastern part of the study area to predict the distribution in the

156

southern part, and vice versa.

we

expected

to

find

a

negative

relationship

between

the

8 157

We conducted an interdisciplinary project to test the above hypotheses, involving

158

sedimentological, hydrodynamic, and ecological expertise. Long-term archived

159

datasets were combined with large, newly recorded datasets of benthos and

160

sediment properties.

161 162

2. Methods

163

2.1 Study area

164

To model the occurrence of Ensis within a coastal gradient, this study was

165

conducted on tidal flats and in sublittoral channels in the eastern Wadden Sea and

166

in the adjacent offshore zone east of 7° 30´ E and south of 55° N, encompassing a

167

total area of about 3,500 km² (Fig. 1). Benthos samples were collected between

168

2009 and 2017. Due to differences in morphology, the study site was subdivided

169

into northern and southern sites (bold rectangles in Fig. 1). The north-eastern part

170

is shaped by islands and barrier islands with deep subtidal creeks with virtually no

171

riverine influence, while in the southern part a few sandy shoals can be found, and

172

the sediment and hydrodynamics are influenced by the estuary of the Elbe river.

173

We therefore anticipated possible differences in the preferred Ensis habitats

174

between these two sites. Furthermore, these two sub-areas were chosen to

175

validate the prediction model for Ensis (see 2.3). We collected 2,393 benthos

176

samples from the northern and 800 samples from the southern study site

177

(subsequently referred to as replicates).

178

9

179 180

Fig. 1: Location of the study area within the eastern Wadden Sea, North Sea (left),

181

and detailed map of the study area (right) subdivided into northern and southern

182

study sites (bold rectangles). Black crosses and red dots indicate the sampling

183

locations and the abundances of American razor clams (E. leei).

184 185

2.2 Sampling design and data recording

186

The abundance and biomass of Ensis were recorded within the whole study area,

187

from the intertidal to the subtidal zone. Intertidal samples were taken using a corer

188

with a sampling area of 107.5 cm², and subtidal samples were taken using a Van

189

Veen grab operated from a research vessel, sampling an area of 1,000 cm². The

190

corer and Van Veen grab both sampled the sediment to a depth of about 20 cm. All

10 191

data recorded with the corer were multiplied by 90.3 to standardise the sampling

192

areas. Samples were rinsed over sieves of 1 mm and 2 mm mesh sizes.

193

Overall, 1,569 stations were sampled, and 3,193 replicates were taken throughout

194

the study area. Three replicates were usually sampled at each sampling station

195

(almost identical geographical coordinates, with slight variation, < 125 m, due to

196

drift of the research vessel). For the modelling approach, we computed the mean

197

values of Ensis abundance and biomass for all replicates per station and used this

198

mean as a dependent variable for the modelling approach (see 2.4).

199

All Ensis individuals collected in the field were stored in plastic bags and kept

200

frozen at −20°C or fixed in 4% formaldehyde for later analysis. All Ensis individuals

201

in each replicate were counted and their length was measured to the nearest mm

202

using calipers. Flesh was separated from the shell and the flesh and shell were

203

weighed separately to the nearest µg. Flesh samples were dried for 24 h at 50°C,

204

the dry weight was determined, and the flesh was then burned in a muffle furnace

205

at 560°C for 12 h. The ash-free dry weight as an indicator of biomass was

206

calculated by subtracting the weight of the ash that remained after burning from the

207

dry weight. The biomass of broken bivalves was determined based on length– or

208

width–mass relationships as previously shown by Dannheim and Rumohr (2012)

209

using our own regressions based on undamaged Ensis individuals.

210 211

2.3 Biotic and abiotic predictors

212

2.3.1 Bathymetry

213

Bathymetric data were provided by the Federal Maritime and Hydrographic Agency

214

at a resolution of 50 m x 50 m for the entire North Sea and German Wadden Sea

215

(Asprion et al., 2013).

11 216 217

2.3.2 Bed shear stress

218

Bed shear stress (τ95; i.e., 95th percentile of local bed shear stress exhibited on 1

219

m² seafloor, expressed in N/m²) was used in the model to estimate the capacity for

220

sediment transport. Current hydrodynamic models provide data of bottom shear

221

stress, which is caused by strong currents and reflects the force of the water on the

222

sediment, providing a direct measure of how benthic organisms might tolerate

223

erosion and re-deposition of the substrate. We evaluated bottom shear stress using

224

a model, which was available on a spatial grid of 1 km x 1 km in the Wadden Sea

225

area and on a grid of 3 km x 3 km in the offshore zone (Kösters and Winter, 2014).

226 227

2.3.3 Bottom shear stress intensity

228

Bottom shear stress intensity (BSSI; i.e., the percentage of time when the bed

229

shear stress exceeds the local critical value for erosion) represents the duration of

230

sediment transport. It is expressed as a dimensionless value ranging from 0 (no

231

time with critical BSSI) to 1 (critical BSSI at all times). As for τ95, BSSI values were

232

taken from Kösters and Winter (2014).

233 234

2.3.4 Duration of inundation

235

The relative time of inundation (in %) during a whole tidal cycle was modelled.

236

Mean relative durations of inundation were calculated for the period from May to

237

June 2011. Model data were obtained by numerical simulations using the coastal

238

ocean model GETM (Burchard and Bolding, 2002; Klingbeil and Burchard, 2013;

239

Klingbeil et al., 2018). The full three-dimensional baroclinic simulations covered the

240

entire Wadden Sea from northern Holland to southern Denmark with a horizontal

12 241

resolution of 200 m x 200 m, 26 adaptive terrain-following vertical layers, and a

242

temporal resolution of 40 s. A detailed description of the model setup and its

243

validation is given in Gräwe et al. (2016). The hydrodynamic data from this model

244

setup were previously used successfully to forecast intertidal seagrass habitats

245

(Folmer et al., 2016). The present study used the simulated sea surface height

246

every 30 min to determine if the tidal flats were inundated or not. Finally, we

247

calculated the mean values for relative duration of inundation for each 200 m x 200

248

m grid point.

249 250

2.3.5 Proportion of mud and fine sand content

251

We collected sediment samples of about 100 g mass in one of the three replicates

252

at each benthos sampling station. We also carried out extensive seafloor

253

recordings (total 1,314 km²) using side scan sonar operated from research vessels

254

in the subtidal area (e.g., Schwarzer and Wittbrodt, 2016; 2018). Additional

255

sediment samples were collected in areas that showed obvious changes in the side

256

scan mosaic. A total of 1,987 sediment samples were collected and the mud (grain

257

size < 63 µm) and fine sand contents (grain size ≥ 63 µm – < 250 µm, Wentworth,

258

1922) were determined in the laboratory. The samples were treated with H2O2 to

259

remove the organic content and dried. Dry weight was recorded, and the samples

260

were rinsed over a cascade of sieves. The sediment remaining in the sieves was

261

dried and weighed and the proportions of the total sediment sample were

262

calculated.

263 264

2.3.6 Co-occurrence with sand mason worm (Lanice conchilega)

13 265

The sand mason worm (L. conchilega) is known to form dense fields in the subtidal

266

area, with up to several thousand individuals per m² (e.g., Rabaut et al., 2009),

267

which can be detected by hydroacoustic devices, especially side scan sonar (e.g.,

268

Schwarzer and Wittbrodt, 2016; 2018). Due to the high density of this polychaete in

269

some areas of the subtidal and its reef-forming nature (Rabaut et al., 2009), we

270

expected it to have negative impacts on the abundance and biomass of Ensis as it

271

has been described before for other benthic species (Ropert and Goulletquer,

272

2000; van Hoey et al., 2008). We therefore recorded the presence and absence of

273

this polychaete in the subtidal zone over a total area of 1,314 km² using side scan

274

sonar (Schwarzer and Wittbrodt, 2016; 2018).

275 276

2.4 Setting up the habitat prediction model

277

Environmental predictors were assigned to a spatial grid of 125 m x 125 m

278

(subsequently referred to as the prediction grid). Bathymetric data were the only

279

predictor available at a higher spatial resolution than the prediction grid. We

280

therefore calculated the mean values of water depth for every grid cell using

281

ArcGIS (version 10.3; Environmental System Research Institute, 2016). For the

282

other predictors, except the proportion of fine sand content, mud and presence of L.

283

conchilega, we used the nearest neighbour value for each of the prediction grid

284

cells. The distribution of fine sand content and mud are based on field data. They

285

were thus only available at a rough spatial scale, and were therefore fitted to the

286

prediction grid with a generalized additive model (GAM; Hastie and Tibshirani,

287

1990) using the function gam() in the R-package 'mgcv' (R Core Team, 2017; R

288

version 3.4.2; Wood, 2017) using latitude and longitude as a smoothed two-

289

dimensional predictor based on cubic splines with the maximal degrees of freedom.

14 290

The result thus represents a cubing interpolation in a given grid cell. This model

291

used a quasibinomial distribution. The predict() function in the 'mgcv' package in R

292

was subsequently used to predict the fine sand values at the coordinates of the

293

prediction grid. The presence or absence of L. conchilega was recorded for all

294

surveyed prediction grid cells. Because this predictor was based on remote sensing

295

(i.e., side scan mosaics; see 2.3.6), no abundance or biomass values were

296

available. As for fine sand content, the probabilities of presence and absence (on a

297

scale of 0–1) were predicted using latitude and longitude.

298

To model the prediction of Ensis, in a first step all predictor values were tested for

299

collinearity using the library usdm (Naimi et al., 2014). Variance inflation factors (vif)

300

were computed. The vif values for water depth and BSSI were >10 and these two

301

predictors were therefore subsequently excluded from the model (Naimi et al.,

302

2014). To set up the Ensis model, we used a GAM with cubic regression splines as

303

smoothing functions and the family attribute “quasipoisson” within the GAM formula

304

(for more details on GAM see Wood, 2017).

305

The following formula was used for the GAM:

306 307

logit (E [Ensisi]) = f1 (τ95i) + f2 (inundation duration) + f3 (proportion fine sandi) + f4

308

(proportion mud) + f5 (presence of L. conchilegai)

(1)

309 310

where Ensis is either the abundance (first model) or the biomass of Ensis (second

311

model), τ95 is the bottom shear stress, inundation duration is the relative time that

312

a prediction grid cell is covered with water, i is the sample number, and f1 to f5 are

313

the smoothing functions. We tested spatial autocorrelation using the function

314

correlog in the R-package “pgirmess” (based on the package “spdep”; Giradoux,

15 315

2018). This revealed no critical spatial autocorrelation for both models and areas

316

(i.e. no autocorrelation for biomass and only a slight autocorrelation for abundance

317

in the distance class of > 3 km for the northern study area which can be

318

disregarded given our small scale modelling grid). Amongst others, the use of GAM

319

is one excellent option to deal with spatial data and habitat predictions. A similar

320

approach to our data, solved with GAMs is presented by Zuur and Ieno (2017).

321

Ensis abundance and biomass were log-transformed to avoid overfitting and to

322

account for some outliers of high values. The model was checked using common

323

selection strategies (i.e. by visual inspection of the residual plots and the function

324

check.gam in the mgcv package; e.g., Korner-Nievergelt et al., 2015; Zuur et al.;

325

2010; 2012). To allow for appropriate ecological interpretation of the results, we set

326

the degrees of freedom to a maximum of k=4 to avoid overfitting (Schwemmer et

327

al., 2016; Schwemmer et al., this issue). Analysis of residuals revealed no violation

328

of linearity, homogeneity, independence, or normality of the random intercept. For

329

each predictor the proportion of explained variance was calculated (e.g. Eskildsen

330

et al., 2013).

331 332

2.5 Prediction of Ensis habitats

333

To validate the Ensis abundance and biomass models, we set up two separate

334

GAMs, one using data from only the northern and the other using data from only

335

the southern sub-area. We compared the performances of both models by

336

predicting the distributions of Ensis abundance and biomass over the entire study

337

area based only on each sub-area, respectively. The prediction was performed

338

using the function ‘predict’ in the R-package mgcv (R Core Team, 2017; R version

339

3.4.2; Wood, 2017). This allowed to predict the habitat suitability for Ensis over the

16 340

entire study area. The raw data for Ensis abundance and biomass were plotted on

341

the interpolated prediction map to visualize the quality of the predictions.

342

Additionally, the spatial distribution of the model residuals was plotted on a map.

343 344

3. Results

345

3.1 Distribution of Ensis

346

Ensis occurred throughout the entire study area, often at high densities (Fig. 1).

347

The maximum density was 80,000 individuals per m², and the maximum biomass

348

was 660.5 g AFDW/m². Ensis densities were highest along a strip off the Wadden

349

Sea islands in the northern study site and in an even broader band in the southern

350

site (Fig. 1). Densities on the intertidal mudflats were much lower, and only

351

increased on low-lying intertidal flats and along the banks of large tidal creek

352

systems in the subtidal or low-intertidal zone (Fig. 1). Virtually no Ensis occurred

353

within the large tidal creek systems themselves. Biomass showed a similar

354

distribution to density except for maximum biomass values in low-lying intertidal

355

flats (particularly in the northern study site), which was mainly attributable to the

356

presence of large adult individuals (see 3.2.4).

357 358

3.2 Ensis habitat model

359

3.2.1 Abundance in the northern study site

360

All five predictors contributed significantly to explaining the abundance of Ensis for

361

the northern study site (Table 1), with an explained variance of 49.2% and a

362

generalised cross validation score (GCV; Wood, 2017) of 0.089. The presence of L.

363

conchilega explained the highest proportion of the variance, while the other

364

predictors provided significant but smaller contributions to the variance (Table 1).

17 365

The relationship between bed shear stress and Ensis abundance in the northern

366

study site followed an optimum curve, with the highest abundances around 1.5

367

N/m² of bed shear stress (Fig. 2a). The probability of encountering Ensis was lower

368

than the predicted mean in areas with lower and higher bed shear stresses. Ensis

369

abundance was low in areas with short periods of inundation and peaked in areas

370

with constant flooding (Fig. 2b). Overall, there was a negative relationship between

371

Ensis abundance and fine sand (Fig. 2c), and an almost linear negative relationship

372

with mud content (Fig. 2d). These findings indicate that Ensis was mainly found in

373

coarser sands with grain sizes ≥ 250 µm. Ensis abundance was positively related

374

to the probability of the presence of L. conchilega up to a probability of 0.6, and

375

then declined rapidly in areas with higher probabilities of L. conchilega occurrence

376

(Fig. 2e).

377 378

Table 1: Modelling results of the GAM for the abundance of Ensis. Relative prop.

379

var. explained = relative proportion of the overall variance explained (49.2% for the

380

northern area; 45.6% for the southern area); edf = estimated degrees of freedom;

381

Ref.df = reference number of degrees of freedom. Northern site Bed shear stress Duration of water coverage Proportion of fine sand Proportion of mud Presence of L. conchilega Southern site Bed shear stress Duration of water coverage Proportion of fine sand Proportion of mud

Relative prop. var. explained

edf

Ref.df

χ²

p-value

2.64 2.87 3.00 2.55 2.95

2.91 2.99 3.00 2.87 3.00

5.54 6.77 3.22 9.86 17.21

0.001 < 0.001 0.021 < 0.001 < 0.001

0.43

edf

Ref.df

χ²

p-value

Relative prop. var. explained

2.893 2.212 2.537 2.916

2.988 2.583 2.86 2.992

8.17 7.648 15.015 4.348

< 0.001 < 0.001 < 0.001 0.007

0.19 0.27 0.41 0.11

0.12 0.22 0.06 0.17

18 Presence of L. conchilega

1.207

1.381

1.552

0.179

382 383

3.2.2 Abundance in the southern study site

384

Four out of five predictors significantly contributed to explaining the abundance of

385

Ensis for the southern study site (Table 1), with an explained variance of 45.6%

386

(GCV=0.037). In contrast to the northern site, L. conchilega had no effect on Ensis

387

abundance (Fig. 2j), while the proportion of fine sand explained the highest

388

proportion of the variance, followed by duration of water coverage, bed shear

389

stress, and the proportion of mud (Table 1).

390

Similar to the northern study site, the relationship between Ensis abundance and

391

bed shear stress in the southern site showed an optimum curve. However, the

392

curve was less pronounced, with a lower maximum (around 1.0 N/m²), and the

393

relationships at the highest bed shear stress values were insecure, as indicated by

394

the broad confidence interval (Fig. 2f). Also similar to the northern site, Ensis

395

abundance was higher in areas with the maximum duration of inundation (Fig. 2g).

396

The maximum Ensis abundance occurred in areas with around 60% of fine sand,

397

followed by a linear decrease at higher values (Fig. 2h). However, in contrast to the

398

northern site, Ensis abundance was high in areas of both low and high mud

399

contents (Fig. 2i).

0.02

19

400 401

Fig. 2 (part 1)

20

402

21 403

Fig. 2: Smoothing functions of GAMs for Ensis abundance in the northern (a-e) and

404

southern study sites (f-j). Smoothed curve: predicted number of Ensis at a given

405

level of the respective predictor variable; dotted red line: overall mean; shaded

406

area: 95% confidence interval; small lines on the x-axis: observation at a given

407

level of the predictor. The smoothed curve and confidence interval above the red

408

dotted line indicates significantly more Ensis than expected according to the

409

predicted mean. If the curve lies below the line, it represents significantly fewer

410

Ensis. Y-axis: relative deviation from the predicted mean of razor clam abundance.

411 412

3.2.3 Biomass in the northern study site

413

All predictors contributed equally to the explained variance for the northern study

414

site, except for the proportion of mud content (Table 2). The overall variance

415

explained was 44.0% (GCV=0.038).

416

Similar to abundance, the biomass of Ensis was optimal in areas with around 1.5

417

N/m² of bed shear stress (Fig. 3a). There was also a steady increase in biomass

418

with increasing duration of inundation (Fig. 3b). As for abundance, biomass was

419

negatively related to the proportion of fine sand, and only increased slightly in areas

420

with much fine sand (Fig. 3c). Ensis biomass was positively related to mud content,

421

but the relationship was not very distinct, as indicated by the broad confidence

422

interval (Fig 3d). Similar to abundance, biomass increased with the probability of

423

encountering L. conchilega, but decreased markedly at high encounter probabilities

424

(Fig. 3e).

425

22 426

Table 2: Modelling results of the GAM for biomass of Ensis. Overall variance

427

explained: 44.0% for the northern area; 41.0% for the southern area. For

428

abbreviations, see Table 1. Northern site Bed shear stress Duration of water coverage Proportion of fine sand Proportion of mud Presence of L. conchilega Southern site Bed shear stress Duration of water coverage Proportion of fine sand Proportion of mud Presence of L. conchilega

edf

Ref.df

χ²

p-value

2.95 1.24 3.00 2.49 2.98

3.00 1.43 3.00 2.83 3.00

52.10 15.38 66.10 5.00 65.60

< 0.001 < 0.001 < 0.001 0.008 < 0.001

Relative prop. var. explained 0.2 0.30

0.25 0.02

0.23

edf

Ref.df

χ²

p-value

Relative prop. var. explained

2.98 1.43 2.47 2.75 2.68

3.00 1.71 2.83 2.95 2.93

15.49 5.25 22.08 15.22 6.70

< 0.001 0.005 < 0.001 < 0.001 < 0.001

0.2 0.24 0.31 0.18 0.07

429 430 431

3.2.4 Biomass in the southern study site

432

The proportion of fine sand explained the highest proportion of the variance for the

433

southern study site, with the other predictors contributing equally, except for the

434

presence of L. conchilega, which had the lowest contribution (Table 2). The overall

435

variance explained was 41.0% (GCV=0.035).

436

As for density, biomass peaked at around 1.0 N/m² of bed shear stress in the

437

southern site (Fig. 3f). Similar to the northern site, biomass increased steadily with

438

increasing duration of water coverage (Fig. 3g). There was a negative relationship

439

between biomass and the proportion of fine sand (Fig. 3h), but Ensis biomass was

440

high in areas with either high and low mud contents (Fig. 3i). Finally, Ensis biomass

23 441

increased slightly with increasing probability of presence of L. conchilega (Fig. 3j),

442

but this relationship was weak, as indicated by the broad confidence interval.

24

443 444

Fig. 3 (part 1)

25

445 446

Fig. 3: Resulting smoothing functions of GAM for Ensis biomass in the northern (a-

447

e) and southern (f-j) study sites. For further explanations, see Fig. 2.

26 448 449

3.3 Predicted Ensis habitats

450

The predicted spatial distributions of Ensis abundance using data from only the

451

northern and only the southern sites, respectively, showed large differences. Using

452

only data from the northern site, the model predicted one hot-spot in the north-

453

western part that matched the raw data well, and a second hot-spot off the

454

Eiderstedt peninsula (for location see Fig. 1) in the central east, which was only

455

partly reflected in the raw data (Fig. 4a). The prediction error of the model was

456

±0.373 individuals / m². In contrast, using data from the southern site predicted a

457

large hot-spot in the northwest (Fig. 4b), and there was a moderately good match

458

between the prediction and the raw data in the south. However, the prediction error

459

for the south model was much higher as compared to the north model (±3.906

460

individuals / m²) and even higher than for a potential overall model (not separating

461

between a northern and a southern site; i.e. 1.213 individuals / m²).

27

462

28 463

Fig. 4: Predicted spatial distribution of Ensis densities throughout the study area,

464

generated using data from the northern (a) and southern study sites (b),

465

respectively.

466 467

Predictions of Ensis biomass distribution using data from the northern site indicated

468

hot-spots in the central north that were well reflected in the raw data, as well as hot-

469

spots in the south, in a line from the Eiderstedt peninsula stretching further south to

470

the Elbe river (Fig. 5a, for location see Fig. 1). The prediction error of the model

471

was ±0.07 g AFDW / m². In contrast, using only data from the southern site

472

predicted a biomass distribution extending further offshore, but with a reasonable

473

match with the raw data both in the vicinity of the Wadden Sea and in the offshore

474

zone (Fig. 5b). The prediction error for the southern site was even less as

475

compared to the north model (±0.036 g AFDW / m²). The prediction error for a

476

potential overall model (not separating between a northern and a southern site)

477

was highest (±0.101 g AFDW / m²).

29

478

30 479

Fig. 5: Predicted spatial distribution of Ensis biomass throughout the study area,

480

generated using raw data from the northern (a) and southern study sites (b),

481

respectively.

482 483

4. Discussion

484

4.1 Distribution of Ensis

485

The distribution of Ensis recorded in this study was in accordance with previous

486

findings of high densities in the (shallow) subtidal zone (Dannheim and Rumohr,

487

2012), in the low-intertidal zone with long periods of inundation (Armonies and

488

Reise, 1999; Dekker and Beukema, 2012), and along the banks of large subtidal

489

creek systems (Houziaux et al., 2012; Gollasch et al., 2015). The maximum density

490

and biomass values recorded in the current study were higher than previously

491

found (Dannheim and Rumohr, 2012; Witbaard et al., 2015), which could be

492

because of the large sample size in the current study or may indicate an increase in

493

the Ensis population.

494

The abundance and biomass of Ensis showed similar distributions, except for a

495

hot-spot in the intertidal zone in the northern study area. This hot-spot represented

496

a low-lying intertidal mud flat with long periods of inundation, which had low

497

numbers (mean: 3.7 individuals / m²) but large biomasses of adults (mean biomass

498

per individual: 1.8 g AFDW) in the corer sample. This was the only location where

499

many adult Ensis were found in the sample. Given that adult Ensis can move up

500

and down their burrows very quickly (e.g., Swennen et al., 1985; Armonies and

501

Reise, 1999), it is possible that the Van Veen grab method might have missed

502

some adult individuals, as indicated by the presence of the upper part of adult

503

Ensis in some grab samples from the subtidal zone. Several samples using a large

31 504

box corer that sampled the sediment to a depth of up to 50 cm revealed whole adult

505

individuals (Fig. 6). The abundance and biomass in the subtidal zone should thus

506

be regarded as minimum values, and the habitat quality for Ensis in this zone is

507

likely to be underestimated. We tried to account for this by also calculating biomass

508

values for broken individuals, but it remains possible that some adult individuals

509

may have been missed in the subtidal zone.

510 511

Fig. 6: Adult Ensis in a sample taken by a box corer in the subtidal (upper picture:

512

zoom on the area indicated with a black box in the lower picture). Several

32 513

individuals could be found at the bottom of their burrows. Note tubes of L.

514

conchilega at the surface.

515 516

4.2 Ensis habitat model

517

Hydrodynamic predictors explained large proportions of the variance in both the

518

abundance and biomass models for Ensis. In accordance with our hypothesis, our

519

results showed that Ensis preferred areas with moderately high bed shear stress,

520

reflecting sediment redistribution. This was in line with experimental studies

521

showing that razor clams may tolerate unstable, coarse, and well-oxygenated

522

sediments (van der Heide et al., 2014), as confirmed by in situ recordings of Ensis

523

distribution patterns (e.g., Dannheim and Rumohr, 2012; Dekker and Beukema,

524

2012; this study). To the best of our knowledge, the current study provides the first

525

statistical evidence to support this relationship over a large area (but see Houizaux

526

et al., 2011). It is likely that the high burying capacity of Ensis (e.g., Swennen et al.,

527

1985; Winter et al., 2012) allows this species to use dynamic sediments with

528

frequent redistributions caused by tidal currents (and maybe also by wave and wind

529

action although these predictors have not been included in the model).

530

Again, in line with our assumptions, our study showed that Ensis preferred areas

531

with long periods of inundation, with peak abundance in the subtidal area. This

532

pattern was confirmed by experimental studies that showed lower growth rates of

533

Ensis on the upper shore line compared with the low-lying intertidal zone, possibly

534

as a result of differences in food supply (Freudendahl et al., 2010). High densities

535

in the low-intertidal and subtidal zones were also confirmed by in situ

536

measurements (e.g., Dannheim and Rumohr, 2012; Dekker and Beukema, 2012;

537

Kottsieper et al., this issue; this study), as well as by ecological niche models

33 538

showing high probabilities of occurrence in water depths up to 30–50 m (Houziaux

539

et al., 2011; Raybaud et al., 2015).

540

As hypothesised, the models for abundance and biomass consistently showed

541

negative relationships with fine sand content for both study sites. However,

542

unexpectedly, the patterns in relation to mud content differed between the study

543

sites, with coarser sediments preferred in the north, but muddy sediments in the

544

south. Virtually all available studies from the North Sea to date have suggested that

545

coarse sediments are the preferred substrate for Ensis (e.g., Dannheim and

546

Rumohr, 2012; Dekker and Beukema 2012). However, preference for muddy sands

547

have been reported at least for parts of the Dutch (Witbaard et al., 2017) and

548

Belgian (Breine et al., 2018) coasts. Additionally, previous laboratory studies

549

showed that Ensis could tolerate high silt contents and may even exhibit higher

550

growth

551

(Kamermans et al., 2013). The current study provides the first evidence for the

552

North Sea across a large area indicating that Ensis can also inhabit areas with very

553

fine sediments. Given that many benthic species are known to prefer higher mud

554

contents (Puls et al., 2012; Singer et al., 2016), these results suggest the potential

555

for competition between Ensis and the native benthic fauna of the Wadden Sea.

556

Interestingly, Ensis was only found in muddy substrates in the south of our study

557

site, particularly in the vicinity of the River Elbe. Although Ensis is known to show

558

lower growth rates at lower salinities (Dannheim and Rumohr, 2012), this may

559

suggest that Ensis might be able to outcompete other benthic species at the lower

560

end of the salinity gradient in the Wadden Sea. Finally, it should be noted that the

561

burying and feeding activities of Ensis itself may increase the amount of mud in the

562

upper layers of the seafloor (Witbaard et al., 2017).

rates

under

these

conditions

compared

with

coarser

substrates

34 563

Our models showed strong effects of L. conchilega on Ensis abundance and

564

biomass, particularly in the northern study area. The probability of encountering

565

Ensis increased with increasing probabilities of encountering L. conchilega up to a

566

threshold, suggesting that the two species shared similar ecological niches.

567

However, the abundance of Ensis declined dramatically when the probability of

568

encountering L. conchilega exceeded 50%, suggesting potential competition. To

569

the best of our knowledge, no studies have examined the direct competition

570

between these two species in the North Sea. However, previous studies found that

571

L. conchilega reduced the carrying capacity for Pacific oysters (Magallana gigas) to

572

19% and was responsible for a 30% oxygen depletion (Ropert and Goulletquer,

573

2000). Despite the different lifestyles of Ensis compared with oysters, they may

574

thus be subject to similar competition. Similar to our findings regarding Ensis

575

biomass and abundance, previous studies indicated that low numbers of L.

576

conchilega were associated with positive effects on species diversity, whereas

577

species richness was decreased in high-density areas (van Hoey et al., 2008).

578

Ecological network analysis (ENA) modelling food web structure and functioning in

579

the study area also suggested competition between these two species (Schückel et

580

al., this issue).

581 582

4.3 Predicted Ensis habitats and quality of the model

583

The prediction models tended to overestimate the abundance and biomass of Ensis

584

throughout the study site (see high number of negative residuals in Fig. 7). While

585

the overestimation was particularly high in areas closer to the shore, the model

586

revealed a slight underestimation of Ensis in the offshore zone. This was true both,

587

for abundance and biomass of Ensis and particularly in the southern study site.

35 588

However, the standard deviations of the residuals were still positive in all cases

589

indicating no significant underestimation of the model.

590

591

36 592

Fig. 7: Spatial distribution of model residuals for (a) abundance and (b) biomass of

593

Ensis. Residuals for an overall model using the entire data set are shown. The

594

residual distribution for the model using only the northern and the southern data

595

set, respectively, are shown in Appendix 1.

596 597

In contrast to our hypothesis, we did not find stable correlations with environmental

598

predictors throughout the entire study area. Data from the northern and southern

599

study sites produced different predicted abundances and biomasses of Ensis. This

600

illustrates the broad ecological niche of this species, and the variation of the

601

realised niche across space. The differences of the realised niche could be

602

attributed to the contrasting morphologies of the northern and southern parts of the

603

Wadden Sea. While hydrodynamic predictors had similar effects at both sites, the

604

presence of L. conchilega only had a strong impact in the north, while the two

605

sediment predictors had the strongest effects in the south. Furthermore, the

606

proportion of fine sand was consistently negatively related to Ensis abundance and

607

biomass at both sites. In contrast, Ensis was found predominately in muddy areas

608

in the south and in coarser substrates in the north. The presence of L. conchilega

609

and sediment composition were thus both likely to be responsible for the different

610

predictions. We therefore concluded that it was necessary to combine the

611

predictions from both sample sites, and the combination of both predictions was

612

shown to produce the greatest overlap between the predicted and observed Ensis

613

abundances and biomasses (Fig. 8). This demonstrates the importance of a

614

modelling approach across a large sea area; the predictions of habitat suitability

615

would have been biased if a prediction for the entire study area had been made

616

based on data from only the northern or only the southern sub-area. The

37 617

development of a universal model for the entire study area was therefore only

618

successful using data from both sampling sites.

619

38 620

Fig. 8: Predicted spatial distribution of Ensis abundance (a) and biomass (b)

621

throughout the study area, generated using raw data from the north for the northern

622

study site and from the south for the southern study site. Both predictions have

623

been combined in GIS.

624 625

4.4 Conclusions

626

In accordance with previous results (e.g., de Mesel et al., 2011; Houziaux et al.,

627

2011; Gollasch et al., 2015), our model showed that Ensis mainly occurred in

628

mobile substrates in the north of our study area (MSRL descriptor D6). Thus, with

629

respect to the descriptors D1 and D4 it can be assumed that competition with other

630

benthic species in these areas is low, given that most other species cannot tolerate

631

these high substrate dynamics (Puls et al., 2012; Singer et al., 2016). However,

632

parts of the Ensis population occurred in muddy sediments, as found in the

633

southern part of our study area (D6). This part of the population could interfere with

634

the native benthic fauna of the Wadden Sea and adjacent offshore zone, which

635

often prefers sediments with high mud contents (Puls et al., 2012; Singer et al.,

636

2016). The degree of overlap with the native fauna is hard to estimate due to a lack

637

of historical data (Houziaux et al., 2011). Although our modelling approach provides

638

indications of a potential overlap of the ecological niche of Ensis with the native

639

fauna (e.g. in case of Lanice), tools like ENA will be needed to disentangle the role

640

of this species in the food web any further. Recent studies suggest that the effect of

641

this neobiotic species on the benthic food web in the intertidal zone is only

642

moderate (Horn et al., 2017; Jung et al., 2017, but see Beukema et al., 2017). In

643

contrast, according to ENA models, there is strong indication of effects on

39 644

biodiversity and food web dynamics in the subtidal (Horn et al., 2017; Schückel et

645

al., this issue).

646

The ecological niche of the Ensis population in the Wadden Sea closely resembles

647

the habitat requirements of this species in its native range (Theroux and Wigley,

648

1983; von Cosel, 2009; Laevitt, 2010). However, in contrast to our study, this

649

species had never been described to occur in muddy sediments within its area of

650

origin. Carnivorous worms and snails are important predators of Ensis in the

651

western Atlantic (von Cosel, 2009). Given that those benthivorous predators are

652

missing in the Wadden Sea and only birds may perform a top-down effect (e.g.,

653

Freudendahl et al., 2010; Tulp et al., 2010; Schwemmer et al., 2016; this issue;

654

Enners et al., 2018), the invasive potential of Ensis in the Wadden Sea (MSRL

655

descriptor D2) is certainly high. Although razor clams in the Wadden Sea grow

656

smaller as compared to their native habitat (von Cosel, 2009), their overall biomass

657

and abundance is much higher. Given the wide spread across many sea areas in

658

Middle Europe a restoration of the system seems impossible (Lackschewitz et al.,

659

2015). If benthivorous predators of Ensis from their native range (von Cosel, 2009)

660

should be introduced to the Wadden Sea, a reduction of the population is possible,

661

however, extinction in the region is unlikely.

662

Finally, morphological differences in the Wadden Sea meant that it was important

663

to perform spatial predictions of abundance and biomass using data from the two

664

sub-areas separately, as the use of just one dataset to predict abundance and

665

biomass for the entire study area would have produced misleading results. Our

666

study thus highlights the need to sample a sufficiently large sea area before

667

predicting the model results to other sites (a prerequisite which is often not possible

668

due to monetary and logistic reasons in ongoing monitoring schemes). However,

40 669

even more importantly, the mismatch of the model results between the northern

670

and southern sub-area is most likely a question of model transferability. Yates et al.

671

(2018) clearly point out that data quality and sampling design are among the most

672

important traits to allow for transfer of models. We tried to resolve these obstacles

673

by collecting a sufficiently large data set using the same sampling design in both of

674

our sub-areas. However, differing morphological conditions in both sites prevented

675

a successful model transfer. Future modelling approaches therefore should aim at

676

finding an appropriate transferability metrics (Yates et al. 2018).

677

Our final model explained a high degree of the variation of the data. However, still

678

more than half of it remains unexplained. One of the most important reasons for

679

this is likely missing additional predictors which also might have hampered model

680

transferability (Yates et al. 2018). We had decided not to increase the number of

681

predictors to avoid overfitting and collinearity (which was already detected testing

682

variance inflation factors among the predictors used). For future modelling

683

approaches, however, it might make sense to include chlorophyll content (as a

684

proxy for the food base of bivalves) as well as salinity (which shows a strong

685

gradient from the Elbe mouth to the northern sub-area and likely might contribute to

686

explaining the differences between the northern and the southern site). Such

687

additional predictors should, however, be available at a similarly small spatial scale

688

as compared to the predictors already used to enable a fine-scale resolution of the

689

prediction. Finally, future studies should aim at combining both, habitat models and

690

ENA. This will allow to link the spatial habitat potential of a species such as Ensis to

691

its effect on the food web at a given site.

692

.

693

41 694

Acknowledgements

695

This study has been carried out in the framework of the projects “STopP-Synthese

696

– from Sediment to Top Predators” (FKZ 03F672B) and “MOSSCO-Synthese –

697

Modular Coupling System for Shelves and Coasts” (FKZ 03F0740B), both funded

698

under the Küstenforschung Nordsee-Ostsee (KüNo) programme of the Forschung

699

für Nachhaltigkeit (FONA) agenda of the German Federal Ministry of Research and

700

Education (BMBF). We are grateful to all the volunteers who helped to collect Ensis

701

samples. We thank the crews of the many research vessels: Ludwig-Prandtl

702

(Helmholtz-Centre Geesthacht), vessels from the Landesbetrieb für Küstenschutz,

703

Nationalpark und Meeresschutz, Egidora, and Littorina from the University of Kiel

704

and Alkor from GEOMAR - Helmholtz Centre for Ocean Research. Finally, we also

705

thank the volunteers and students who helped with the laboratory work. S. Furness

706

provided linguistic support.

707 708

Declaration of competing interests

709

The authors declare that they have no competing interests.

710 711

42 712

References

713

Armonies, W., Reise, K., 1999. On the population development of the introduced

714

razor clam Ensis americanus near the island of Sylt (North Sea). Helgoländer

715

Meeresuntersuchungen 52, 291-300.

716

Asprion, U., Sbresny, J., Griffel, G., Elbracht, J., 2013. Die Bathymetrie der

717

deutschen Nordsee–Erstellung der projektweiten Bezugsfläche. Report of the

718

Geo-scientific Potential of the German North Sea (GPDN) Project. Federal

719

Institute for Geosciences and Natural Resources (BGR), the State Authority for

720

Mining, Energy and Geology of Lower Saxony (LBEG) and the Federal

721

Maritime and Hydrographic Agency (BSH). https://www.gpdn.de/?pgId=404.

722

(last accessed 22. January 2019).

723

Baird, D., Asmus, H., Asmus, R., 2012. Effect of invasive species on the structure

724

and function of the Sylt- Rømø Bight ecosystem, northern Wadden Sea, over

725

three time periods. Mar Ecol Prog Ser 462: 143-161.

726

Beukema, J.J., Dekker, R., Drent, J., van der Meer, J., 2017. Long-term changes in

727

annual growth of bivalves in the Wadden Sea: influences of temperature, food,

728

and abundance. Marine Ecology Progress Series 573, 143–156.

729

Burchard, H., Bolding, K., 2002. GETM - a General Estuarine Transport Model.

730

Scientific Documentation. Technical Report EUR 20253 EN. European

731

Commission.

732 733 734 735

Buschbaum, C., Lackschewitz, D., Reise, K., 2012 Nonnative macrobenthos in the Wadden Sea ecosystem. Ocean & Coastal Management 68, 89-101. Breine, N.T., de Backer, A., van Colen, C., Moens, T., Hostens, K., van Hoey, G., 2018.

Structural

and

functional

diversity

of

soft-bottom

microbenthic

43 736

communities in the southern North Sea. Estuarine, Coastal and Shelf Science

737

214, 173-184.

738 739 740 741

Cadée, G. C. 2000. Herring gulls feeding on a recent invader in the Wadden Sea, Ensis directus. Geological Society, Special Publication 177, 459–464. Dannheim, J., Rumohr, H., 2012. The fate of an immigrant: Ensis directus in the eastern German Bight. Helgand Marine Research 66: 307–317.

742

Dekker, R., Beukema, J.J., 2012. Long-term dynamics and productivity of a

743

successful invader: the first three decades of the bivalve Ensis directus in the

744

western Wadden Sea. Journal of Sea Research 71, 31-40.

745

Enners, L., Schwemmer, P, Corman, A.-M., Voigt, C.C., Garthe, S. 2018.

746

Intercolony variations in movement patterns and foraging behaviors among

747

herring gulls (Larus argentatus) breeding in the eastern Wadden Sea. Ecology

748

and Evolution https://doi.org/10.1002/ece3.4167

749 750

Environmental Systems Research Institute (ESRI), 2016. ArcGIS v. 10.3 ESRI. Redlands, California.

751

Eskildsen, A., le Roux, P.C., Heikkinen, R.K., Høye, T.T., Kissling, W.D., Pöyry, J.,

752

Wisz, M.S., Luoto, M., 2013. Testing species distribution models across space

753

and time: high latitude butterflies and recent warming. Global Ecology and

754

Biogeography 22, 1293–1303.

755

Feldens, P., Schulze, I., Papenmeier, S., Schönke M., Schneider von Deimling, J.,

756

2017. Improved Interpretation of marine sedimentary environments using multi-

757

frequency

758

doi:10.3390/geosciences8060214

multibeam

backscatter

data.

Geosciences

2018,

8,

214;

44 759

Folmer, E.O., van Beusekom, J.E.E., Dolch, T., Gräwe, U., van Katwijk, M.M.,

760

Kolbe, K., Philippart, C.J.M., 2016. Consensus forecasting of intertidal seagrass

761

habitat in the Wadden Sea. Journal of Applied Ecology 53, 1800-1813.

762

Freudendahl, A.S.L., Nielsen, M.N., Jensen, T., Jensen, K.T., 2010. The introduced

763

clam Ensis americanus in the Wadden Sea: field experiment on impact of bird

764

predation and tidal level on survival and growth. Helgoland Marine Research

765

64, 93-100.

766

Galil, B.S., Marchini, A., Occhipinti-Ambrogi, A., Minchin, D., Narščius, A., Ojaveer,

767

H., Olenin, S., 2014. International arrivals: widespread bioinvasions in

768

European Seas. Ethology Ecology & Evolution 26, 152-171.

769

Giraudoux, P., 2018. pgirmess: Spatial Analysis and Data Mining for Field

770

Ecologists.

771

project.org/package=pgirmess (last accessed 22. January 2019)

772 773

R

package

version

1.6.9.

https://CRAN.R-

Gollasch, S., Nehring, S., 2006. National checklist for aquatic species in Germany. Aquatic invasions 1: 245-269.

774

Gollasch, S., Kerckhof, F., Craeymeersch, J., Goulletquer, P., Jensen, K., Jelmert,

775

A., Minchin, D., 2015. Alien Species Alert: Ensis directus. Current Status of

776

Invasions by the Marine Bivalve Ensis directus. ICES Cooperative Research

777

Report 323: 1–32.

778

Gräwe, U., Flöser, G., Gerkema, T., Duran-Matute, M., Badewien, T.H., Schulz, E.,

779

Burchard, H., 2016. A numerical model for the entire Wadden Sea: Skill

780

assessment and analysis of hydrodynamics. Journal of Geophysical Research -

781

Oceans 121, 5231-5251.

782 783

Hastie, T., Tibshirani, R., 1990. Generalized Additive Models. Chapman and Hall, BocaRaton, Florida.

45 784

Houziaux, J.S., Craeymeersch, J., Merckx, B., Kerckhof, F., van Lancker, V.,

785

Courtens, W., Stienen, E., Perdon, J., Goudswaard, P.C., van Hoey, G., Vigin,

786

L., Hostens, K., Vincx, M., Degraer, S., 2011. 'EnSIS' - Ecosystem Sensitivity to

787

Invasive Species. Final Report. Brussels: Belgian Science Policy Office 2012 –

788

Research Programme Science for a Sustainable Development. 105 pp.

789

http://www.belspo.be/belspo/ssd/science/Reports/ENSIS%20FinalReport_%20

790

ML.pdf (last accessed 14th February 2019).

791

Horn, S., de la Vega, C., Asmus, R., Schwemmer, P., Enners, L., Garthe, S.,

792

Binder, K., Asmus, H., (2017) Interaction between birds and macrofauna within

793

food webs of six intertidal habitats of the Wadden Sea. PloOne DOI:

794

10.1371/journal.pone.0176381.

795

Jung, S., Dekker, R., Germain, M., Philippart, C.J.M., Witte, J.I.J., van der Veer

796

H.W., 2017. Long-term shifts in intertidal predator and prey communities in the

797

Wadden Sea and consequences for food requirements and supply. Marine

798

Ecology Progress Series 579, 37–53.

799

Kamermans, P., Brummelhuis, E., Dedert, M., 2013. Effect of algae- and silt

800

concentration on clearance- and growth rate of the razor clam Ensis directus,

801

Conrad. Journal of Experimental Marine Biology and Ecology 446, 102–109.

802

Katsanevakis, S., Inger, W., Zenetos, A., Leppäkoski, E., Çinar, M.E., Oztürk, B.,

803

Grabowski, M., Golani, D., Cardoso, C., 2014. Impacts of invasive alien marine

804

species on ecosystem services and biodiversity: a pan-European review.

805

Aquatic Invasions 9, 391-423.

806

Klingbeil, K., Burchard, H., 2013. Implementation of a direct nonhydrostatic

807

pressure gradient discretisation into a layered ocean model. Ocean Modelling

808

65, 64-77.

46 809

Klingbeil, K., Lemarié, F., Debreu, L., Burchard, H., 2018. The numerics of

810

hydrostatic structured-grid coastal ocean models: state of the art and future

811

perspectives. Ocean Modelling 125, 80-105.

812

Korner-Nievergelt, F., Roth, T., von Felten, S., Guelat, J., Almasi, B., Korner-

813

Nievergelt, P., 2015. Bayesian data analysis in ecology using linear models

814

with R, BUGS, and Stan. Academic Press, 328 pp.

815 816

Kösters, F., Winter, C., 2014. Exploring German Bight coastal morphodynamics based on modelled bed shear stress. Geo-Mar Letters 34, 21-36.

817

Lackschweitz, D., Reise, K., Buschbaum, C., Karez, R., 2015. Neobiota in

818

deutschen Küstengewässern. Eingeschleppte und kryptogene Tier- und

819

Pflanzenarten an der deutschen Nord- und Ostseeküste. Landesamt für

820

Landwirtschaft, Umwelt und ländliche Räume des Landes Schleswig-Holstein

821

(LLUR),

Gewässer;

D25.

822

https://www.umweltdaten.landsh.de/nuis/wafis/kueste/neobiota.pdf

(last

823

accessed 17th June 2019).

LLLUR

SH

–

824

Leavitt, D.F., 2010. Biology of the Atlantic jacknife (razor) clam (Ensis directus

825

Conrad, 1843). Aquaculture Centre, University of Maryland. NRAC Publication

826

No.

827

https://pdfs.semanticscholar.org/458b/39aead771289253fa46cafa332c961a7ed

828

34.pdf?_ga=2.208759888.372541613.1560759788-2025400810.1560759788

829

(last accessed 17th June 2019).

217-2010.

830

Naimi, B., Hamm, N.A.S., Groen, T.A., Skidmore, A.K., Toxopeus, A.G., 2014.

831

Where is positional uncertainty a problem for species distribution modelling.

832

Ecography, 37, 191-203.

47 833

Puls, W., van Bernem, K.-H., Eppel, D., Kapitza, H., Pleskachevsky, A.,

834

Riethmüller, R., Vaessen, B., 2012. Prediction of benthic community structure

835

from environmental variables in a soft-sediment tidal basin (North Sea).

836

Helgoland Marine Research 66, 345-361.

837

R Core Team, 2017. R: A language and environment for statistical computing. R

838

Foundation

for

Statistical

Computing,

Vienna,

839

project.org/. (last accessed 22th January 2019).

Austria.

https://www.R-

840

Rabaut, M., Vincx, M., Degraer, S. 2009. Do Lanice conchilega (sandmason)

841

aggregations classify as reefs? Quantifying habitat modifying effects. Helgoland

842

Marine Research 63, 37–46.

843

Raybaud, V., Beaugrand, G., Dewarumez J.-M., Luczak, C. 2015. Climate-induced

844

range shifts of the American jackknife clam Ensis directus in Europe. Biological

845

Invasions 17, 725–741.

846 847 848 849

Reise, K., 1985. Tidal flat ecology - an experimental approach to species interactions. Springer, Berlin. Reise, K., Olenin, S., Thieltges, D.W., 2006. Are aliens threatening European aquatic coastal ecosystems? Helgoland Marine Research 60: 77-83.

850

Ropert, M., Goulletquer, P., 2000. Comparative physiological energetics of two

851

suspension feeders: polychaete annelid Lanice conchilega (Pallas 1766) and

852

Pacific cupped oyster Crassostrea gigas (Thunberg 1795). Aquaculture 181,

853

171-189.

854

Sahlin, U., Ryden, T., Nyberg, C.D., Smith, H.G., 2011. A benefit analysis of

855

screening for invasive species – base-rate uncertainty and the value of

856

information. Methods in Ecology and Evolution 2, 500-508.

48 857

Schwemmer, P., Güpner, F., Adler, S., Klingbeil, K., Garthe, S., 2016. Modelling

858

small-scale foraging habitat use in breeding Eurasian oystercatchers

859

(Haematopus ostralegus) in relation to prey distribution and environmental

860

predictors. Ecology Modelling 320, 322-333.

861

Schwarzer, K., Wittbrodt, K., 2016. ALKOR – 440 – Cruise report, 30.06.-

862

12.07.2014. STopP I (Vom Sediment zum Topp-Prädator – Einfluss von

863

Eigenschaften des Meeresbodens auf Benthos und benthivore Vögel,

864

Teilprojekt StopP-See). DOI: 10.3289/CR_AL_440 (last accessed 14th February

865

2019).

866

Schwarzer, K., Wittbrodt, K., 2018. LITTORINA 17-14 North Sea – Dithmarschen

867

Bay. Cruise report. KÜNO-Projekt STopP. DOI: 10.13140/RG.2.2.31143.80808

868

(last accessed 14th February 2019).

869 870

Severijns, N. 2002. Distribution of the American Jack-knife clam Ensis directus in Europe 23 years after its introduction. Gloria Maris, 40, 63–111.

871

Singer, A., Schückel, U., Beck, M., Bleich, O., Brumsack, H.-J., Freund, H.,

872

Geimecke, C., Lettmann, K.A., Millat, G., Staneva, J., Vanselow, A., Westphal,

873

H., Wolff, J.-O., Wurpts, A., Kröncke, I., 2016. Small-scale benthos distribution

874

modelling in a North Sea tidal basin in response to climatic and environmental

875

changes (1970s-2009). Marine Ecology Progress Series 551, 13-30.

876

Swennen, C., Leopold, M.F., Stock, M., 1985. Notes on growth and behavior of the

877

American razor clam Ensis directus in the Wadden Sea and the predation on it

878

by birds. Helgoländer Meeresuntersuchungen 39, 255-261.

879

Theroux, R.B., Wigley, R.L., 1983. Mollusks based on specimens in the National

880

Marine Fisheries Service Woods Hole Collection. NOAA Technical Report

49 881

NMFS

SSRF-768.

https://spo.nmfs.noaa.gov/sites/default/files/legacy-

882

pdfs/SSRF768.pdf (last accessed 17th June 2019).

883

Tulp, I., Craeymeersch, J., Leopold, M.F., van Damme, C., Fey, F., Verdaat, H.,

884

2010. The role of the invasive bivalve Ensis directus as food source for fish and

885

birds in the Dutch coastal zone. Estuarine, Coastal and Shelf Science 90: 116-

886

128.

887

Van der Heide, T., Tielens, E., van der Zee, E.M., Weerman, E.J., Holthuijsen, S.,

888

Eriksson, B.K., Piersma, T., van der Koppel, J., Olff, H. 2014. Predation and

889

habitat modification synergistically interact to control bivalve recruitment on

890

intertidal mudflats. Biological Conservation 172, 163-169.

891

Van Hoey, G., Guilini. K., Rabaut, M., Vincx, M., Degraer, S., 2008. Ecological

892

implications of the presence of the tube-building polychaete Lanice conchilega

893

on soft-bottom benthic ecosystems. Marine Biology 154, 1009-1019.

894 895 896 897 898 899

Von Cosel, R., 2009. The razor shells of the eastern Atlantic, part 2. Pharidae II: the genus Ensis Schumacher, 1817 (Bivalvia, Solenoidea). Basteria 73: 9-56. Wentworth, C.K., 1922. A scale of grade and class terms for clastic sediments. Journal of Geology. 30: 377-392. Winter, A.G., Deits, R.L.H., Hosoi, A.E., 2012. Localized fluidization burrowing mechanics of Ensis directus. Journal of Experimental Biology 215, 2072-2080.

900

Witbaard, R., Duineveld, G.C.A., Bergman, M.J.N., Witte, H.I.J., Groot, L.,

901

Rozemeijer, M.C.J., 2015. The growth and dynamics of Ensis directus in the

902

near-shore Dutch coastal zone of the North Sea. Journal of Sea Research 95,

903

95–105.

50 904

Witbaard, R., Bergman, M.J.N., van Weerlee, E., Duineveld, G.C.A., 2017. An

905

estimation of the effects of Ensis directus on the transport and burial of silt in

906

the near-shore Dutch coastal zone of the North Sea. 127, 95–104.

907 908

Wood, S.N., 2017. Generalized Additive Models: An Introduction with R. London: Chapman and Hall, 2nd edition, 416 pp.

909

Yates, K.L., Bouchet, P.J., Caley, M.J., Mengersen, K., Randin, C.F., Parnell, S.,

910

Fielding, A.H., Bamford, A.J., Ban, S., Barbosa, A.M., Dormann, C.F., Elith, J.,

911

Embling, C.B., Ervin G.N., Fisher, R., Gould, S., Graf, R.F., Gregr, E.J., Halpin,

912

P.N., Heikkinen, R.K., Heinänen, S., Jones, A.R., Krishnakumar, P.K. and 27

913

others, 2018. Outstanding challenges in the transferability of ecological models.

914

Trends in Ecology and Evolution 33: 790-802.

915 916 917 918

Zuur, A.F., Leno, E.N., Elphick, C.S., 2010. A protocol for data exploration to avoid common statistical problems. Methods in Ecology and Evolution 1, 3–14. Zuur, A.F., 2012. A beginner’s guide to Generalized Additive Models with R. Highland Statistics Ltd, Newburgh, 206 pp.

919

Zuur, A.F., 2017. Beginner´s guide to spatial, temporal and spatial-temporal

920

ecological data analyses with R-INLA – Volume II: GAM and zero-inflated

921

models. Highland Statistics Ltd, Newburgh, 363-708 pp.

Spatial distribution of model residuals for abundance of Ensis using two separate models for the northern and the southern study site.

Spatial distribution of model residuals for biomass of Ensis using two separate models for the northern and the southern study site.

Highlights: • • • • •

Long-term data were used to determine ecological niche and predict habitat potential of Ensis leei Abundance and biomass of Ensis were related with sedimentologic and hydrodanamic predictors Ensis occurred both, in coarse sediments as well as in muddy environments Ensis abundance and biomass declined in areas with dense Lanice conchilega reefs The invasive Ensis is thought to overlap with the pristine benthic fauna particularly in muddy sediments