A risk-based surveillance design for the marine pest Mediterranean fanworm Sabella spallanzanii ( Gmelin , 1791 ) ( Polychaeta : Sabellidae ) – a New Zealand case study

To determine the presence, or extent and spread, of marine pests is often difficult and decisions on allocating limited sampling effort need to be made using available information. This study presents a robust structured methodology to develop a detection survey for the marine pest Sabella spallanzanii. The design of the detection survey used modelled hydrodynamics of the area and expert knowledge on settlement characteristics for Sabella. Habitat suitability for settlement was defined based on expert opinion elicited using the Analytical Hierarchy Process (AHP) technique and a self-administrated questionnaire. Zones for Sabella settlement were then identified by overlaying suitable habitat areas and hydrodynamic patterns of potential larval propagule dispersal. Settlement zones were assigned a risk/likelihood ranking to ensure available surveying effort was allocated efficiently over a potentially wide settlement area. This design was shown to be successful in detecting Sabella. Provided underlying hydrodynamic information is available, the structured approach to pest species detection presented here could readily be applied to develop surveillance plans for other broadcast spawning marine pests.


Introduction
The Mediterranean fanworm (Sabella spallanzanii (Gmelin, 1791)), hereafter referred to as Sabella, is recognised as an invasive marine organism that can have negative impacts where it becomes established (Holloway and Keough 2002;O'Brien et al. 2006;Murray and Keable 2013;Ross et al. 2013).Extensive literature exists describing these impacts and a few examples are given here.Sabella alters the structure of the habitat by forming a dense canopy of filter feeding structures above the long tubes in which the animal is housed, thus causing competitive interactions in understory taxa.The dense canopy and efficient filter feeding can also cause changes to the water flow rates which has effects on food availability and feeding efficacy of other organisms and inhibiting their recruitment (Holloway and Keough 2002).Impacts in aquaculture are known with the potential for this species to become a nuisance fouler, where a cultured species may need to compete with Sabella for nutrients and space (Murray and Keable 2013).
The Mediterranean fan worm, as the name suggests, is native to the Mediterranean and Atlantic coast of Europe (Currie et al. 2000).It has now been recorded from Brazil where it is thought to be an introduction (Currie et al. 2000), Australia (Murray and Keable 2013) and the first occurrence in New Zealand was in 2008 (Read et al. 2011).Sabella is an unwanted and notifiable organism under the New Zealand Biosecurity Act 1993.In 2008 a single large, mature individual of Sabella was collected from a wharf pile in Lyttleton harbour during a diver survey (Figure 1; Read et al. 2011).Subsequently, a second, larger population of Sabella was detected in Auckland's Waitemata harbour, after which eradication attempts were abandoned but surveillance for spread continued (Figure 1; Read et al. 2011).
This fouling species is now well established in this northern region of New Zealand where commercial and recreational vessels are likely to act as vectors spreading this species to other regions across New Zealand.In 2013, Sabella was first reported in the Coromandel region, on two infested vessels (Figures 1  and 2).At this time the organism was only recorded as being on these vessels, not as being established on the substrate in the area.This instigated a joint agency, local and central government response with an initial delimitation survey, an in-water hull cleaning process and a 6-month follow up detection/delimitation survey in case of accidental propagule release during vessel cleaning.
After action was taken to clean the vessels of all the mature tubeworms the vessels were relocated back to the Auckland region where Sabella is well established.Although the cleaning method followed strict protocols developed in the field (unpublished) to prevent potential release of Sabella into the environment, there is a chance that some adults or larvae were released into the harbour.A follow up survey was undertaken six months post-defouling to ascertain whether there were any Sabella that may have spawned and/or "escaped" and had established on the surrounding substrate.Six months was considered the appropriate amount of time, as any new recruits that had settled would then be of a detectable size (i.e., detectable by naked eye of a diver) (Murray and Keable 2013).If Sabella were detected, decisions needed to be made about whether to continue to try to eradicate or control this population.
This paper describes the design of this follow up survey to maximise the chance of detection of Sabella recruits using a structured, pragmatic and targeted approach incorporating local conditions as well as expert opinion.This design was enacted and was successful in detecting Sabella.

Methods
To develop a robust surveillance plan with the objective of detecting Sabella if it is present, with a high level of confidence, we decided to use stratified random sampling within a defined area.This approach proportionately allocated more sampling effort to the "high risk" strata.These strata consisted of habitats that were more likely to favour Sabella settlement.

Sampling Extent defined by hydrodynamics
In order to define a realistic sampling frame, an outer boundary was determined by modelling the potential natural larval spread from the release sites (the vessel locations) under similar environmental conditions experienced during the cleaning activities.Specific meteorological conditions were observed as calm through to storm conditions during the modelled period.In particular spawning conditions included a period where "… the wind and swell direction seemed to be predominantly from the NW with winds in excess of 30 knots for a sustained period of time over a 24-hour period." (pers. comm. D. Hodges WRC 26 July 2013).
A calibrated 3-D hydrodynamic model for the Firth of Thames and Hauraki Gulf (Knight and Beamsley 2013) was used to simulate this spread scenario.The model was based on the Semi-implicit Eulerian-Lagrangian finite-element model (SELFE; Zhang and Baptista 2008) framework.The SELFE model is an unstructured model, composed of a range of triangular-mesh sizes from 20 m close to shore to 4.5 km in offshore areas.The model has been shown to perform well at reproducing tidal elevations, water temperature structure and velocities at five sites around the Firth of Thames over a 6-week period (Knight and Beamsley 2013).
Larval dispersal was modelled using abiotic neutrally-buoyant Eularian tracers using the default advection transport scheme in the SELFE model.This approach allowed for the continuous release of "larvae" for a month from both the surface and bottom cells at the two vessel locations.
A month-long period of virtual larvae releases was modelled beginning 23 February 2009, with mean concentration statistics calculated from the resulting tracers calculated for a 7-day period, beginning two weeks into the model run (i.e.9-20 March 2009).The modelled period covered the range of tidal and wind conditions that occurred during the time the vessels were in the Coromandel Harbour.The twoweek period was selected to approximate the pelagic larval duration of the species (14 to 21 days; NIMPIS 2013), with two storm events occurring near the beginning of the modelled period.
In order to estimate the modelled tracer concentrations in propagule concentration units, tracer concentrations were scaled based on an estimated spawning population size on the vessels and the number of propagules per individual spawning Sabella.
The spawning population size was based on estimates of the surface area of each barge, the density of adult animals and an estimated proportion of the population that would be able to spawn on each vessel.The surface area of each vessel was estimated at 274 m 2 .This was calculated based on a horizontal dimension of 25 m × 6 m and subsurface depth of 2 m for each barge.The density of Sabella was assumed to be 1000/m 2 on the most heavily fouled vessel and 800/m 2 on the other vessel.The proportion of spawning adult Sabella organisms was assumed to be 40% of the total population.Consequently the total number of spawning individuals was estimated to be 109 600 and 87 680 on each vessel.
The number of propagules was estimated at 50 000 per spawner per spawning event.Two spawning events were assumed to occur during a 2-week period, which implied a spawning rate of 0.0827 propagules per second per spawner (G.Read pers.comm).
When this information was combined with the spawning population estimates, this gave propagule supply rate estimates for the 2-week period of: • 9 064 propagules/second from vessel one • 7 251 propagules/second from vessel two This information was used to scale arbitrary tracer concentrations from the model and ensure the units presented in the results could be interpreted as a propagule concentration (i.e.propagules/m 3 ).
The scaled model outputs of both the mean and maximum propagule concentrations for the weeklong period were passed to the ArcGIS 10.1 software for spatial analysis.The concentration results were then segregated into four categories using an unbiased Jenks "natural breaks" classifier in the ArcGIS 10.1 software.
The Jenks optimisation approach creates categories that reduce variance within groups, while maximising variance between groups.This has the effect of partitioning the data into classes based on natural groups, or breaks, based on low points in the data distribution.We chose the boundary line created by these natural breaks that fitted with a spatial extent that was practicable to survey within real world constraints.Hence, the final sampling extent was defined by using the boundary of larval propagule concentration in the hydrodynamic dispersion model, plus the addition of a few high settlement likelihood habitats (three marine farms and one reef) that fell just outside of this boundary.The entire intertidal area was excluded as previous studies (CSIRO 2001;NIMPIS 2013) suggested Sabella was restricted to subtidal environments.

Defining high settlement likelihood habitats by Expert Elicitation using an Analytical Hierarchy Process
An expert elicitation exercise was used to establish which regions in the Coromandel harbour could be considered the high risk strata.As mentioned previously, "high risk" or "high likelihood" is defined as habitat on which Sabella is often found established and known to settle on.As with most marine invasive species, there is limited understanding of the ecological and environmental requirements and tolerance of Sabella in newly occupied/invaded habitats (Holloway and Keough 2002;O'Brien et al. 2006).In order to improve our understanding of habitat suitability for this species, an expert elicitation exercise was conducted with 14 experts (including both New Zealand and international experts).All participating experts had experience with sabellids and in marine biosecurity.Eight experts had specific experience with Sabella spallanzanii.
The elicitation exercise was designed using the Analytical Hierarchy Process (AHP, Saaty 1987) as is described in Correa-Berger (2007), Kendrick and Saaty (2007) and Saaty (2008).A questionnaire using pairwise comparisons where experts would rank the importance of the environmental variables (in this case depth, sediment type and exposure) was distributed to the panel via email.Experts also needed to use pairwise comparisons to rank the categories considered within each variable.Experts' pairwise answers were combined using the geometric mean and integrated in an AHP matrix (Saaty 1987(Saaty , 2008)).The final (i.e.combined expert) weight, or importance, for each variable was also estimated using the geometric mean.Based on the rankings the various habitat types received from this exercise, a spatial analysis was performed where these habitat layers were overlaid with the sampling area defined by the hydrodynamic dispersion model to identify zones where Sabella would be most likely to settle.

Sample size calculations
To establish the appropriate number of transects to be searched to maximise chances of detecting any Sabella, should it be present, methodology as described in Smith (2006) was used.Smith's (2006) method calculates the probability of detecting a "rare" species based on search area, assumed search efficiency and assumed density.It has been demonstrated to work well at a range of spatial distributions and low densities.
To calculate the area that must be searched one uses the below formula: Where β is the estimated search efficiency and µ is the estimated species density.
Solve for search area ∝ at the appropriate confidence level ( C ).
As per the methods described in Smith (2006), we tried a range of confidence levels, search efficiencies and assumed maximum densities of Sabella as a "rare" occurrence.

Sample units
The sampling unit recommended was a 100 m 2 transect.Transects are a commonly used method for diver visual searches (McCormick and Choat 1987;Kingsford 1998;Willis et al. 2000;Pande and Gardner 2009) and this is also a practicable size when breaking down a survey area into smaller searchable units.

Proportional allocation of sampling units to high settlement likelihood habitat
A grid of 100 m × 1 m cells was generated in ArcGIS which were designed to simulate a 100 m long diver transect survey with a visibility width of 1 meter.This diver survey grid was overlaid on the defined sampling area.Those survey cells that were not suitable for dive surveys (e.g.located on land or in the intertidal areas) were excluded from the sampling frame.Each included grid point was then assigned a priority value, according to which habitat classification the centre point of the grid intersected.Even if one transect spans several different types of habitat, the priority ranking will remain constant if the centre points of all grids are used to assign this weight.
The "Sampling Priority Value" (SPV) was defined as: The weights (W) were estimated by the AHP and assigned to each variable as described previously.
Each cell was then assigned a priority value.The weighting is equivalent to a probability that Sabella will be found in that cell.Accordingly the cells can be grouped and the appropriate number of transects randomly selected from each priority group.In order to select the cells from each priority group for the surveillance a code in the R statistical software package (R Core Team 2015), which factors in the weighting (priority group) of each cell, was used.The "R" code applies these probabilities sequentially so that choosing the next item is proportional to the weights amongst the remaining items.This ensures that more sampling effort goes into the higher probability (priority value) cells.

Sampling extent defined by Hydrodynamics
The results of the model show that the highest mean and maximum propagule concentrations both occur within Coromandel Harbour and concentrations are reduced by at least a factor of 10 outside of the harbour (Figure 3).Therefore we used the mean propagule concentrations over the last week of the month-long simulation to define the sampling area (Knight 2013).We chose the Jenks' classification break (1.57-3.12log 10 propagules/m 3 ) which fitted closely to the natural harbour boundaries and covered an area that was practicable to survey.

Expert Elicitation (AHP)
The results of the AHP exercise showed that experts considered the type of sediment, followed by exposure, as the most important determining variable for Sabella (Figure 4).Breaking this down more, the AHP results showed that of the "sediments", experts considered reefs and man-made structures (e.g., farms, moorings) localised in shallow (<20 m) sheltered environments as high settlement likelihood area for this species.

Sample size
Using Smith's (2006) equation as described previously, various values of search efficiency and assumed density were trialled.The range of options considered are presented in Table 1.Assuming a 0.2 search efficiency (probability of detecting an individual given it is within the search area) where the assumed conservatively low density was 0.001/m 2 at a confidence level of 95% (95% probability of detecting it if it is present) -the total search area is 14 978 m 2 .
This was considered to be an area that was practical to search given resource constraints.This sampling area translated into 150 sampling units (100 m 2 transects).These transects were then proportionally allocated to the high risk/high settlement likelihood habitats as per Figure 5.

Discussion
The methodology presented here describes a robust way of developing a detection survey for a marine pest in an area where it was previously believed absent.Surveillance in the marine space is complex, as it is difficult to access this environment and there are no true boundaries within which one can easily define a surveillance area (Corsin et al. 2009;Pande et al. 2015).Hydrodynamic modelling to determine  This figure shows the ranking assigned to each of the variables in the Analytical Hierarchy process show sediment type to be the most important variable (A) , of which man-made structures, followed by biogenic reef and reef to be the most important (C).For the other variables sheltered areas (B) and shallow depths (D) were also considered more important for the settlement of Sabella.
Table 1.Different options for transect numbers to search the Coromandel harbour for Sabella.The highlighted row is the one that was considered the most appropriate to employ for this survey taking into account practically and constraints.This is assuming a 0.2 search efficiency with an assumed density of 0.001/m 2 at a confidence level of 95%, yielding a result of a 14978 m 2 area to search.the dispersal range of a species of interest is one way to help define boundaries in the marine environment (Pande et al. 2015).Finding a robust method for pest detection is of paramount importance for decision makers who need to identify whether it is feasible or not to try and eradicate the pest in question, or whether some other control method might be more appropriate.The methodology used in this study was an effort to maximise the probability of detection, within real world constraints.Data for the variable "exposure", calculated as mean current velocity at a site, were not found to be an adequate representation of exposure, hence were not considered in the analysis.However, the sensitivity to detect Sabella is not meaningfully lowered by excluding these data, as the AHP showed the predilection for Sabella to establish on man-made structures to be the most important factor, with depth also playing an important role.Exposure may be of assistance for estimating suitability if more appropriate indicators of exposure (e.g. from wave models) were available.
The usefulness of AHP in decision making for environmental impact assessments and capturing and analysing socio-economic impacts has been demonstrated (Ramanathan 2001).More recently the use of AHP, or expert opinion has been seen in the biosecurity space (Simon et al. 2016;Soliman and Inglis 2016;Burgman 2016).AHP is particularly useful when decisions need to be made using sparse or incomplete information.It involves consulting with a range of experts and standardises aggregation of their judgements to avoid bias (Ramanathan 2001).Since decisions about eradication or management of an invasive species often need to be made when little is known about the true potential impact of the organism in a new environment, expert opinion almost always comes into play.AHP is a structured way of utilising expert opinion and it also allows the relative importance of various factors to be considered when ranking the options (Kendrick and Saaty 2007).The main criticism of AHP has been "rank reversal" where the addition of another option can change the relative position of the original choices; however when the concept of absolute measurement is used ("high", "medium", "low") for the rankings, rank reversal does not occur (Ramanathan 2001).This was the approach used in the AHP this study.
The survey design suggested in this study shares some fundamental aspects with the CRIMP (Centre for Research on Introduced Marine Pests) baseline port survey methodology (Hewitt and Martin 2001).The CRIMP methodology focusses on sampling habitats and sites that are the most likely to be colonised by the target species and that any local knowledge about where it might be found should be taken into consideration.Hewitt and Martin (2001) also operate under the assumption that they are most likely to be detected near the point of inoculation.However, this study has gone further through the application of an AHP, hydrodynamic modelling and proportional allocation of sampling effort to high settlement likelihood sites to objectively apply allocate sampling effort.Hewitt and Martin (2001) aimed to detect rare species at a density of 0.1 individuals at 95% proba-bility, and to achieve this, suggested sampling 13 units, or a minimum of 7 units.In this case the surveyed units were berths in a port environment.Given the Coromandel harbour is not a busy port that the CRIMP methodology was designed for, it was appropriate that we used different sample units (i.e.diver survey transects) and assumed a lower density of the target species, resulting in a higher number of sample units.
This robust surveillance design was applied to increase the likelihood of detection, and even though the total area surveyed in proportion to the total potentially infested area is small, Sabella individuals were detected using this methodology.There is always a possibility of missing pest recruits, simply because of the extremely large area to be surveyed.However, in most cases where surveillance for a pest is undertaken, similar constraints are likely to apply.However, the aim of this surveillance design was detection and this goal was achieved when the surveillance was enacted.
Other, more simple and targeted niche surveillance activities were also enacted at the same time as the survey described here, to allow for detection of Sabella individuals that may have come from sources other than larval recruits originating from the two infested vessels in question.Thus the structured surveillance approach described here became one of a suite of multiple different surveillance methods, all of which together provided a comprehensive way to detect Sabella in this area.

Conclusions
The potential area of infestation of pest incursions in marine environments, are often much larger than can be practicably surveyed by available resources.Consequently, we conclude that structured approaches to surveying, which allow for targeted sampling effort, and based on the best available knowledge, are a necessity.

Figure 1 .
Figure 1.Map of New Zealand showing the original infestation sites of Sabella in Lyttleton Harbour, Waitemata Harbour in Auckland and the location of this study in Coromandel Harbour.

Figure 2 .
Figure 2. Coromandel harbour bathymetry, expressed as metres below mean sea level (MSL) and barge release locations (marked with "+") that were used in the spread modelling exercise.

.
: Weight of Category X within Variable i : Weight Variable i n: Total number of variables considered

Figure 3 .
Figure 3. Mean (left) and maximum (right) tracer concentration results for a week-long period occurring after 14-18 days of continuous tracer release.The colours correspond to tracer concentrations displayed on a logarithmic scale.

Figure 4 .
Figure 4.This figure shows the ranking assigned to each of the variables in the Analytical Hierarchy process show sediment type to be the most important variable (A) , of which man-made structures, followed by biogenic reef and reef to be the most important (C).For the other variables sheltered areas (B) and shallow depths (D) were also considered more important for the settlement of Sabella.

Figure 5 .
Figure 5. Shows the boundaries of the sampling area overlaid with the high risk habitats and excluding the intertidal habitat (light green).The boundary is defined by the Jenks' natural break category.The high risk areas are as the highest weighted risk categories being darker colours.The sampling points are proportionally allocated to the high risk categories.