Distribution of the invasive Caprella mutica Schurin , 1935 and native Caprella linearis ( Linnaeus , 1767 ) on artificial hard substrates in the North Sea : separation by habitat

Studying offshore natural and artificial hard substrates in the southern North Sea (51oN–57oN/1oW–9oE), the invasive introduced Japanese skeleton shrimp Caprella mutica Schurin, 1935 was found to co-exist with the native Caprella linearis (Linnaeus, 1767) only on near-shore locations that had an intertidal zone (e.g., wind farm foundations). In contrast, on far offshore and strictly subtidal locations, such as shipwrecks and rocky reefs, only C. linearis was found. Based on these exploratory observations, we hypothesised that artificial structures that are only subtidal are inhabited exclusively by C. linearis, and never by C. mutica. To test this hypothesis and understand factors driving each species’ habitat preferences, habitat suitability models were constructed using generalised additive models, based on samples collected in 2013–2015 from offshore gas platforms, buoys, shipwrecks, and rocky reefs and combined with data from other published and unpublished surveys (2001–2014). The models showed that the presence of C. mutica is explained by the availability of intertidal and floating hard substrates, suspended particulate matter density (SPM), mean annual sea surface temperature, salinity, and current velocity. The C. linearis model included subtidal hard substrates, SPM, salinity, temperature, and current velocity. The modelled distributions showed a significant difference, demonstrating that C. linearis’ habitat preference does not fully overlap with that of C. mutica. Thus, the native and alien Caprella species are likely to be able to co-exist in the North Sea.


Introduction
Non-indigenous invasive species can threaten native species and cause extinctions (Simberloff 2010).In the marine environment, no extinctions have been directly linked to invasions, but invasive species may still drive native species to potentially less preferred habitats (Hill and Lodge 1999).This distribution shift can prevent extinction of the native species (Gurevitch and Padilla 2004) but alters the local species community.
Fast-growing introduced species are likely to establish viable populations after introduction (Sakai et al. 2001).The invasive Japanese skeleton shrimp Caprella mutica Schurin, 1935, a caprellid amphipod, has fast reproductive capabilities (Cook et al. 2007b;Shucksmith et al. 2009).Caprella mutica females can release their first brood within 24-26 days after they themselves hatch, allowing a rapid expansion once introduced (Cook et al. 2007b).In its native range in north-east Asia (Schurin 1935;Ashton 2006), it is found associated with macro-algae in shallow water (Fedotov 1991;Ashton et al. 2007) and on aquaculture structures (Kawashima et al. 1999).Caprella mutica was first recorded in Europe in the Netherlands in 1994, and only 14 years later it had been identified on all coasts along the North Sea, English Channel, and Celtic Sea (Platvoet et al. 1995;Cook et al. 2007a).Using mitochondrial DNA, Ashton et al. (2008) noted that C. mutica may have been initially introduced to European waters in Scotland and Norway, with secondary dispersal from there.Such secondary dispersal is most likely aided by flotsam and human activities, e.g., ship traffic and aquaculture (Ashton 2006).This amphipod is known to cling to ship hulls, which may be an important factor for longer range transportation (Cook et al. 2007a;Frey et al. 2009;Adams et al. 2014;Zabin et al. 2014).Caprella mutica is able to survive 20 days without food, enabling high survival during long transportation periods (Cook et al. 2007b).
Both C. mutica and C. linearis are observed on anthropogenic structures (Ashelby 2005;Buschbaum and Gutow 2005;Page et al. 2006;Cook et al. 2007a;Bouma and Lengkeek 2013;Macleod 2013;Nall et al. 2015).In the North Sea, large numbers of artificial structures are present, all providing potential habitat for both species.Several authors have suggested that offshore anthropogenic structures can function as stepping stones for invasive species (Mineur et al. 2012;Adams et al. 2014;De Mesel et al. 2015).
Stepping stones provide habitat in an environment normally unsuitable for survival (MacArthur et al. 1967), facilitating species to spread faster to potential recipient locations than they otherwise could.Adams et al. (2014) suggests C. mutica may be transported by ships between the coast and offshore structures, which they then colonise.From there, the invasive species may spread even further.Shipwrecks (strictly subtidal) are common artificial structures in the North Sea (Coolen et al. 2015b).Many other man-made structures are also present in the North Sea, including buoys, wind farms and offshore oil and gas (O&G) installations.One difference between most shipwrecks and the other structures is that the others penetrate the surface, providing surface / inter-tidal habitats.All these artificial structures represent potential habitat for C. mutica and C. linearis.In previous studies of fouling diversity on artificial structures in the North Sea, only C. linearis was recorded from shipwrecks (Zintzen and Massin 2010;Lengkeek et al. 2013) while both species were recorded from turbine foundations of wind farms in Dutch waters (Bouma and Lengkeek 2013).This suggested that C. mutica may be excluded from artificial structures that are strictly subtidal whereas these could be colonised by C. linearis.To test this hypothesis, we reviewed available sampling data and conducted additional sampling programs of fouling communities in the North Sea.With advances in remote sensing and modelling techniques (Gayer et al. 2006;Brown et al. 2011), it was possible to obtain datasets with sufficiently fine resolution to model the combined effects of a range of variables on the presence of species (Reiss et al. 2014).In this study, we investigated the distribution of C. mutica and C. linearis in the offshore southern and central North Sea using habitat suitability models based on generalised additive models (GAM) to determine if the habitats occupied differed.

Study area
The North Sea (surface area of 575,300 km 2 ) is a coastal sea in the north east Atlantic Ocean.It is largely surrounded by land with tidal water entering from the Atlantic Ocean via the English Channel to the south and between the United Kingdom and Norway to the north (southern and central North Sea: Figure 1).Water circulation in the North Sea is driven counter-clockwise by tides and wind.Sea surface temperatures vary between 2-8 °C in winter and 12-21 °C in coastal waters during summer (Otto et al. 1990;Rijkswaterstaat 2015).The input of nutrients from rivers in the surrounding countries has a strong influence on turbidity in coastal zones (Brockmann et al. 1990; Figure 2).
The seafloor of the North Sea largely consists of sandy and muddy (mobile) sediments, interrupted by zones of coarse material, mainly gravel and rocks (Coolen et al. 2015a).The area covered with coarse

Acquisition of samples
We collected 150 samples from 64 unique locations in the Dutch and British Exclusive Economic Zones of the North Sea.Samples were extracted from different sources, between 2013 and 2015.
-Mixed macrofauna samples were obtained from offshore gas production installations, shipwrecks, the Borkum Reef Grounds, and Texel Rough.We used surface supplied diving equipment (platforms) and scuba equipment (all other objects).We collected samples using an airlift sampler fed by surface supplied air or by air from side-mounted scuba tanks.A 500 cm 2 frame was attached to the sampled surface and the surface scraped clean using a putty knife.The suction part of the airlift was held close to the putty knife to collect the scraped specimens, which were deposited in a 500 µm mesh size net.The methods used are described in more detail in Coolen et al. (2015a); -In addition to airlift sampling, we actively searched for caprellids on additional shipwrecks during wreck inventories in 2013-2015.During these dives we searched for caprellids for at least 15 minutes, or until at least 25 specimens were collected from different surfaces.Caprellids are easily spotted under water due to their elongated bodies protruding from the surfaces they cling to.
We followed methods also applied in Coolen et al. (2015b).About 50 dives (including airlift sampling) were made in water depths between 20 and 45 meters between 2013 and 2015; -During the diving expeditions described above, we did not survey every location visited, but inspected ghost fishing nets removed from the ship wrecks by other scuba divers for presence of caprellids at some locations; -During inspection, repair and maintenance (IRM) work by SSE divers on offshore gas installations, caprellids clinging to divers suits were collected opportunistically from the suits after the IRM divers submerged.Both species tend to cling to clothes and working gear when removed from the substrate (pers.obs.J.W.P. Coolen); -During IRM work on offshore buoys in Dutch waters, qualitative macrofauna samples were taken from the marine growth left on deck after cleaning of the buoys.At all buoy locations, about 1 litre of mixed macrofauna with mussels was collected by the IRM vessel crew and frozen in a household freezer; and -Eight macrofauna samples were obtained from a gravel reef and the surrounding sandy sediments in the Borkum Reef Grounds using a box-corer (0.076 m 2 ).Only samples with >15 cm penetration depth were retained.Detailed methods are described in Coolen et al. (2015a).
All samples were stored in a borax-buffered formaldehyde-seawater solution (6%) except the buoy samples, which were frozen in a −20 °C freezer.All observed caprellids were identified to the lowest taxonomic level possible, using the World Register of Marine Species (WoRMS Editorial Board 2015) as a standard for taxonomical nomenclature.For identification, keys by Stock (1955), Larsen (1998) and Guerra-García (2014) were consulted.

Additional presence-absence data
Further data with presence-absence observations of C. mutica and C. linearis were acquired from published and unpublished sources.In order to ensure that only true absence data were included, confirmation that both species were sought was asked from the people involved in creating the data.If such confirmation was not possible, the species was excluded from the model data for that specific dataset (NA in supplementary Table S1).If the identification of both species was unsure, the full dataset was excluded.The following sources were included in the analysis: - A dataset of all locations is provided as online supplement and the data obtained from newly sampled locations were deposited in the Dryad Digital Repository (Coolen et al. 2016).

Model creation
The predictor variables used to model the presence and absence of C. mutica and C. linearis (Table 1) were obtained from various sources and included mean annual sea surface temperature (MASST), salinity (MASSS), current velocity (MASSCV), suspended particulate matter (SPM) and three types of hard substrates; subtidal (SHS), intertidal (IHS) and floating hard substrate (FHS).The categorical hard substrates location data were converted to presence (1) and absence (0) raster data (presence-absence, PA).Combination into a single factor variable was not possible since offshore installations placed on the sea bottom contain both subtidal and intertidal substrates.Sea bottom depth was initially considered as a possible predictor variable but was not included in the models because the available data only included bottom depth and did not reflect the true local depth, e.g. on intertidal or floating locations.Data for all species were combined into a single dataset (Table S1).Within this dataset, all observations were converted to PA data to compensate for differences in sampling methods.Using ArcGIS 10.2.1.3497for Desktop (ESRI, Redlands, CA), both species presence-absence and predictor variable data were fitted to 1 km square grid cells extending between 50ºN and 61ºN, and 4ºW and 9ºE.A subset of unique locations from an area outside the dataset was excluded at this stage but used for later validation of the optimised models.This subset comprised the data from the buoy and wind farm surveys in Belgium and the September 2015 wreck survey.Every grid cell containing an observation of presence or absence of either C. mutica or C. linearis was exported, including the value of each predictor variable per cell.This resulted in a training dataset with PA observations of 160 locations for C. mutica and of 137 locations for C. linearis.
For analysis, R version 3.2.0(R Core Team 2015) and RStudio version 0.98.1103(RStudio 2014) were used.Data were explored following methods described in Zuur et al. (2010).Using Cleveland dotplots, boxplots, pairplots, Pearson correlation coefficients, variance inflation factors, and multi-panel scatterplots (Cleveland 1985;Sarkar 2008;Dormann et al. 2013) the presence of outliers, multi-collinearity, relationships and interactions were analysed for combinations of all caprellid observations and predictor variables.Species distribution often shows a non-linear relation with continuous abiotic data such as temperature (Austin 2007).Therefore, a GAM was created for both species, using the gam function in the mgcv package (Wood 2011).The data were modelled using the Bernoulli distribution with logit link function.GAMs are prone to overfitting (Wood and Augustin 2002); therefore, the number of knots for the smoothers for the continuous variables was limited to four.All continuous variables were initially included as smoothers and the need for smoothing evaluated during model optimisation.The initial full model took the following form [s(*) indicates smoothed variables]:

~
To exclude predictor variables and optimise the models, backward selection using Akaike Information Criteria (AIC; Akaike 1973) was combined with ecological evaluation of the influence of included and excluded effects.The optimal models were validated to verify the underlying assumptions.Model residuals were plotted against fitted values to analyse homogeneity of variance and against all covariates used during model selection to assess model fit.The part of the data that was excluded at the beginning of the analysis contained 69 locations for C. mutica and 63 for C. linearis.This dataset was used to validate the predictions of the models using a generalised linear model (GLM).The relation between C. mutica and C. linearis was then modelled using a GLM and C. mutica presence was added to the validated C. linearis model to evaluate this relation in combination with the predictor values, following Ros et al. (2015).
With the resulting models and the full dataset of predictor variables, a presence absence prediction raster was calculated using the predict.gamfunction (Wood 2011).To test for differences between the distributions of C. mutica and C. linearis, the presence absence predictions of both species were compared using a Students T-Test for paired samples.The resulting presence absence prediction grids for both species were visualised using ArcGIS.

Obtained data
From 289 possible locations C. mutica was present 74 times and C. linearis 41 times.Within these locations, C. mutica and C. linearis only co-existed in two wind farm locations near the Dutch coast (Table S1; Figure 1).Caprella mutica was only present on floating objects or intertidal offshore structures.C. mutica was never observed on any of the subtidal shipwrecks or rocky reef locations.Most C. linearis observations were from offshore locations with intertidal or subtidal hard substrates with a few observations on floating hard substrates.On most wrecks, C. linearis was present as well as on most O&G installations, and the rocky reefs of the Borkum Reef grounds, and the Texel Rough.Caprella mutica and C. linearis only co-occurred at intertidal locations.Neither species was observed on any of the sandy sediment locations.

Model selection and visualisation
During model selection for mutica, all predictor variables from the initial full model were included in the final model.All the continuous data (temperature, salinity, current velocity and SPM) were left in the model as non-linear terms (Figure 4).The resulting model to predict the presence and absence of C. mutica in the North Sea was: and subtidal hard substrates (p = 0.059) were nonsignificant but had explanatory value and were therefore kept in the model.The deviance explained by the model was 57.2% with an adjusted r 2 of 0.517.During model selection for C. linearis, floating hard substrates were removed since inclusion of this variable resulted in a strong negative effect, and AIC evaluation forced subsequent removal of subtidal hard substrates from the model.The removal of floating hard substrates (keeping subtidal hard substrates) made more ecological sense.After removal of floating hard substrates, the AIC score for the model with subtidal hard substrates was lower than without.Intertidal hard substrates was removed as this had a low explanatory value resulting in a lower AIC for the model without intertidal hard substrates.The smoothers for both salinity and SPM were linear (edf = 1).Temperature and current velocity were included as non-linear terms (Figure 5).The resulting model to predict the presence and absence of C. linearis was: . ~ 66.6783 1.3725 * 1.9492 * 0.4642 * Significant effects were found for subtidal hard substrates (p = 0.036), temperature (p <0.001), salinity (p = 0.008) and SPM (p = 0.020).Contrary to C. mutica, SPM had a negative effect on the presence of C. linearis.Current velocity was nonsignificant (p = 0.231) but did explain part of the deviance and so was kept in the model.The deviance explained by the model was 49.9% with an adjusted r 2 of 0.529.
The predicted distributions differed significantly (t-test, t = 27.875, p <0.001).The relation between C. mutica and C. linearis was non-significant (GLM, p = 0.529) and inclusion of C. mutica presence in the C. linearis model resulted in a higher AIC than the validated model without C. mutica as well as a nonsignificant effect for C. mutica (p = 0.111).
The training set for the C. mutica model included data from the northern North Sea but the C. linearis set did not.Furthermore, no other data were obtained from the North Sea area north of the Dogger Bank.Therefore the model validation and visualisation were limited to the southern and central North Sea (51ºN-57ºN / 1ºW-9ºE).Validation showed the model for C. mutica explained 37.5% of the deviance of the validation dataset (GLM, p = 0.007) and for C. linearis it explained 20.3% of the deviance (GLM, p = 0.005).This is 66% and 40% of the deviance explained by the main models, respectively.The model visualisations showed C. mutica was highly suited to nearshore waters and partially suited to artificial structures that were floating or had intertidal surfaces in more offshore waters (Figure 6).Caprella linearis was highly suited to both nearshore to offshore waters and was suited to locations with completely subtidal reefs (Figure 7).

Discussion
The modelled distributions of the alien, invasive, C. mutica and native C. linearis showed a significant difference.Caprella mutica occurred more frequently nearshore whereas C. linearis occurred more frequently offshore.Suspended particulate matter (SPM) density, which includes detritus, the main food source for both species (Guerra-García and Tierno de Figueroa 2009), had a contrasting effect on the species.SPM had a significant positive effect on the presence of C. mutica and a negative effect on the presence of C. linearis.This suggests that C. linearis prefers habitats with lower levels of SPM, such as waters that are far offshore (Brockmann et al. 1990;Jickells 1998).
There was a clear negative non-linear relation between current velocity and presence of C. mutica.Macleod (2013) also demonstrated the negative impact of high currents on C. mutica.Current velocity was included in the model for C. linearis, but there was no significant effect.This is in line with findings by Macleod (2013) who attributes this to the smaller body size of C. linearis, mean adult length 5.12 mm, compared to C. mutica, mean adult length 11.39 mm (Shucksmith et al. 2009).However, intertidal sites are more exposed to other hydrodynamic factors such as wave induced currents, resulting in higher disturbance than is present at subtidal sites (England et al. 2008).Guerra-García (2001) observed that caprellid species on wave-exposed sites were larger than those on sheltered sites.These findings contradict with our results and should be explored in future models, by including, for example, mean wave agitation as an additional factor (e.g., Dutertre et al. 2013).
Both caprellid species showed a significant nonlinear relation with mean annual sea surface temperature, with an optimum between 11 and 12 °C for C. linearis and an optimum at lower temperatures, <11 °C, for C. mutica.Boos (2009) reported a significant increase of C. mutica mortality when experimentally kept at continuous temperatures >20°C.Caprella mutica may have a preference for the lower range of average temperatures in the North Sea which may explain the limited reports of C. mutica presence from southern Europe (Cook et al. 2007a but see Almón et al. 2014), where mean annual sea surface temperatures are between 19 and 21 °C, with over 25 °C for summer averages (Skliris et al. 2012).Caprella linearis is also known from northern waters with lower temperatures than observed here (Larsen 1998).Shucksmith et al. (2009) suggest temperature as well as salinity play a role in the co-existence of both species.Salinity had a small but significant effect on C. linearis, while C. mutica showed a non-linear relation that may partly be explained by the lack of data from the central to northern North Sea, where the salinity is between 34.5 and 35.The model visualisation (Figure 6) was limited to the area with salinity <34.5.
Notable was C. mutica's association with shallow water objects, which meant they had a high potential to encounter rafts, such as macroalgae.Both C. mutica and C. linearis are known to use such rafts (Thiel and Gutow 2005).Rafting may have aided in the rapid dispersal of C. mutica in European waters.Buschbaum and Gutow (2005) propose that C. mutica may have colonised Helgoland using rafts, and Ashton  (2006) showed this species' ability to use drifting algae for dispersal over distances > 5 km.Caprella linearis may have colonised intertidal locations in a similar manner.However, to colonise strictly subtidal reefs, rafting does not help since drifting objects are unlikely to reach a sub-tidal-only location.The colonisation of these locations may have been aided by objects transported by currents near the bottom, such as free rolling sponges that have been observed in the southern North Sea (pers.obs.J.W.P. Coolen).A combined preference for shallow waters and high SPM density, explains why C. mutica was absent from intertidal hard substrate that was far offshore while C. linearis was present there.In the current study C. mutica was never observed at depths >17 m.Caprella mutica is known from depths between 0 and 17 m and C. linearis between 0 and 65 m (Fedotov 1991;Moen and Svensen 2004;Vanagt and Faasse 2014).Our observations confirm that C. linearis is able to occupy deeper locations, such as wrecks, which are outside C. mutica's preferred depth range.Offshore hard substrates are scarce in the North Sea.Therefore, artificial hard substrates likely play an important role in the distribution of C. linearis in offshore waters.
Conflicting with the absence of C. mutica from deep subtidal hard substrates was that our models showed that C. mutica presence was partly explained by the presence of subtidal hard substrates.This resulted because near-shore installations were modelled as being both intertidal and subtidal habitats.Our models had a high "predictor to response variable" ratio, which may cause overfitting (Babyak 2004) and decrease model accuracy (Vittinghoff and Mcculloch 2006;Wisz et al. 2008).However, the models were validated with an independent dataset that demonstrated their ability to predict the presence of the modelled species.The deviance explained by these predictions was limited to 37.5% for C. mutica and 20.3% for C. linearis.To improve these percentages, more data are needed and additional explanatory variables such as mean wave agitation need to be explored.Further model improvement could include using not only mean annual values for predictor variables but also annual maxima or minima, for example for sea surface temperatures.We limited the prediction area to the southern and central North Sea because only a limited amount of data were available for the northern parts.Extrapolations to other areas should be validated using additional survey data.
Habitat selection is dependent on competition between species and competition between species is habitat dependent (Morris 2003): strong competitors in one habitat may be less competitive in another habitat.This was shown by Ros et al. (2015) for the relationship between other native and non-native caprellids around the Iberian Peninsula.There, the non-native species was dominant in more southern waters that were warmer and saltier.In the present study, no data on C. linearis distribution prior to and immediately after the introduction of C. mutica to the North Sea were available.We were therefore unable to draw conclusions about if competition from C. mutica had displaced C. linearis from former habitats.In any case, hard deep water habitats, including artificial structures, represent refuges for C. linearis, where populations can survive away from possible competition from C. mutica.

Figure 1 .
Figure 1.Caprellid distributions in the southern and central North Sea.Locations with non-native Caprella mutica only: red dots, with native Caprella linearis only: yellow dots, with both species: black dots, with neither species: white dots.

Figure 3 .
Figure 3. Map with locations of all artificial structures (black dots) used for modelling.Depth contours as in Figure 1.
Macrofauna surveys at wind farms: Dutch -Egmond aan Zee (Bouma andLengkeek 2013) and Princess Amalia Wind Farm(Vanagt and Faasse 2014), Danish -Horns Rev 1(Leonhard and Frederiksen 2006), German -Alpha Ventus and FINO 1 research platform(Krone et al. 2013;Gutow et al. 2014) and Belgian -C-Power and Belwind (De Mesel et al. 2015); -Invasive species surveys around the Shetland Islands (unpublished data supplied by Shucksmith and Shelmerdine, NAFC Marine Centre, 2015); -Surveys of wrecks in Belgian waters (Zintzen and Massin 2010); -Post drill surveys near the German A6-A installation on the Dogger Bank (Glorius et al. 2014a) and the German L1-2 installation in the Borkum Reef Ground area (Glorius et al. 2012); -Baseline survey for a proposed pipeline between the German A6-A installation and Danish Ravn platform (Glorius et al. 2014b) and a proposed exploration well in the German B11-5 block on the Dogger Bank (Glorius et al. 2013); and -Buoy surveys in Belgian waters between 1998 and 2014 (Kerckhof, unpublished data).

Figure 4 .
Figure 4. Predicted relations between environmental variables and probability of Caprella mutica presence for: mean annual sea surface temperature (MASST; °C), mean annual sea surface salinity (MASSS), mean annual sea surface current velocity (MASSCV; m.s -1 ), and mean annual density of suspended particulate matter (SPM; g.m -3 ).Dashed lines show 95% confidence interval, short vertical lines on the x axis show the density of samples in the model.

Figure 5 .
Figure 5. Predicted relations between environmental variables and probability of Caprella linearis presence for: mean annual sea surface temperature (MASST; °C), and mean annual sea surface current velocity (MASSCV; m.s -1 ).Dashed lines show 95% confidence interval, short vertical lines on the x axis show the density of samples in the model.
Within this model, significant effects were found for temperature (p = 0.043), salinity (p = 0.015), current velocity (p = 0.014) and SPM (p = 0.013), with SPM showing a positive effect in most of the data range.Effects of floating (p = 0.059, intertidal (p = 0.141)

Table 1 .
List of variables used in the habitat suitability models for Caprella mutica and Caprella linearis with description of the variable, units used, study period, and data source.