Uncovered pathways: Modelling dispersal dynamics of ship-mediated marine introduced species

1. Marine traffic is the main vector for marine non-indigenous species (NIS) that may hitchhike in ballast water tanks or attached to vessel hulls. Understanding ma rine traffic dynamics and estimating the associated risk of NIS dispersal reveals points of leverage for preventive NIS management. This study presents a method to quantify the risk of ship hull fouling-mediated NIS dispersal and identifies main dispersal hubs in marine traffic networks. We use the Galapagos Marine Reserve (GMR) as a case study to test the applicability of this method. 2. Ship position data derived from the automatic identification system (AIS) served as a basis for a network consisting of nodes (moorings and anchorages) and edges (ship routes). Wetted surface areas (WSA) describe the parts of vessel hulls


| INTRODUC TI ON
Marine traffic is the main anthropogenic pathway for marine non-indigenous species (NIS; Carlton, 1979;Molnar et al., 2008;Shucksmith & Shelmerdine, 2015). Once established in a new suitable habitat, marine NIS may have negative effects on local ecosystems by outcompeting native species, and may ultimately lead to a decrease in biodiversity and cause economic losses (Bax et al., 2003;Pyšek & Richardson, 2010). Vessels transport species as stowaways in their ballast water tanks or attached to their hulls (Godwin, 2003;Ruiz et al., 1997). With an ever-increasing world fleet-mainly due to the growing marine cargo sector (UNCTAD, 2018)-the risk of NIS dispersal is growing. In addition to large commercial vessels, recreational and passenger vessels represent an important vector group, dispersing marine NIS at a more local dimension (Ashton et al., 2014;Ferrario et al., 2017;Iacarella et al., 2020;Ulman et al., 2019).
Eradication measures are very costly in time and money, and often fail to eliminate the unwanted marine NIS (Sambrook et al., 2014). Therefore, management should intervene at an early stage, focusing on preventing the arrival of marine NIS and developing quick response strategies (Hulme, 2009;Ojaveer et al., 2015). With more available data about vessel movements and the use of modern modelling tools, science is able to provide policymakers with crucial information to develop such regulations (Laeseke et al., 2020).
Modelling the risk of NIS introduction associated with hull fouling has been addressed by previous scientific studies. In one of the few examples, Floerl et al. (2009) modelled the spread of marine NIS as a consequence of vessel traffic between marinas in New Zealand, and found that busy marinas are more likely to become infected by NIS. However, both less frequented and busy marinas contributed to the spread of NIS. Another example, based on ship voyages between a major port and provinces in Indonesia, identified those with the highest numbers of arriving ships and hence exposed to the largest risk of marine NIS introductions (Azmi et al., 2015). Finally, the most recent example applied graph theory to detect most intensely used shipping routes between NIS hotspots and marine-protected areas (MPA), identifying recreational vessels as the most important vessel group contributing to the connectedness between them (Iacarella et al., 2020). In all mentioned examples, the models did not consider varying vessel sizes, meaning that vessels contributed equally to NIS dispersal. However, several studies have shown that vessel dispersal capacities depend on individual vessel characteristics (Floerl & Inglis, 2005;Verling et al., 2005). Moser et al. (2016) developed a method to determine the 'wetted surface area' (WSA) of a vessel, which is the area of the hull underneath the water line. In other words, the WSA describes the area available for fouling species to settle on the ship hull and therefore may be used as a proxy for individual vessel dispersal capacity.
The Galapagos Marine Reserve (GMR) is a hotspot for recreational and touristic vessel traffic, as the annual visitor numbers surpassed 270,000 in 2018 (OTG, 2020). Current biosecurity regulations prohibit the discharge of ballast water in the GMR and focus on vessel fouling of entering ships. Yet, the number of known marine NIS in Galapagos has recently increased from 6 (Keith et al., 2016) to 53, as a result of a thorough analysis of fouling assemblages (Calder et al., 2019;Carlton et al., 2019;Keppel et al., 2019;Lambert, 2019;McCann, 2019).
Possible reasons might be a recent increase in artificial structures nearby Galapagos' ports  and a growth in touristic yachting in Galapagos in the last decades (OTG, 2020).
In this study, we combine a marine traffic network based on individual vessels' WSAs and tracks with a numeric model to quantify dispersal capacities. Changes of overall dispersal capacity caused by sequential removals of network elements (nodes, edges and ship classes) indicated their importance. We use the GMR as a case study to identify hotspots and key linkages of NIS dispersal and translate our results into preventive NIS management recommendations.
Focusing solely on the pathway component, this method is applicable to any place with available ship position data.

| Study area
The Galapagos Islands, Ecuador, are located in the Pacific Ocean, on the equator, about 1,000 km west of the Pacific coast of Latin America ( Figure 1). Five of the islands are inhabited; Santa Cruz (15,701), San Cristobal (7,199), Isabela (2,344), Floreana (143;INEC, 2010) and Baltra (Ecuadorian Navy and Airforce bases).
For this study, we adapt the terminology of the ports so that they reflect the name of the island: Puerto Ayora (Port Santa Cruz),

| Ship data
The Automatic Identification System (AIS) was originally developed for ship collision avoidance, and is now required for all commercial vessels above 300 gross tonnes on international voyages and all passenger vessels by the International Maritime Organization (IMO).
The AIS is composed of transceivers installed on vessels, as well as satellite and terrestrial receivers. The Galapagos National Park provided AIS data covering the GMR in 2016. The data were recorded at least every few seconds and involved the following relevant K E Y W O R D S automatic identification system, dispersal modelling, Galapagos, invasive species, network analysis, non-indigenous species, risk management variables: ship identification number, time stamp, latitude, longitude and speed. The raw data were cleaned and processed in several steps to suit our analysis using the R software (R Core Team, 2019) and the geosphere package (Hijmans, 2019; see Appendix S1 for details). We added information about vessels, such as length, tonnage and ship type from the homepage www.marin etraf fic.com using web scraping in R. We grouped all vessels in the following ship classes: Recreational (private yachts), passenger (day-tour boats and liveaboard cruisers), fishing, general cargo, tanker, container, other (patrol, research and other cargo) and unknown. We identified anchorages or moorings (from here on called anchorages) for all vessels based on their speed and geographical position, and matched those with marine visitor sites and ports in the GMR (DPNG, 2014). All vessels entering or leaving the GMR were assigned to the node Outside, a placeholder node located north-east of the reserve.

| Wetted surface area
We calculated the WSAs of vessels following the method outlined in Moser et al. (2016) using either ships' tonnage, length or width, depending on their availability. For vessels with a volume smaller than 115 gross tonnes, the method of Moser et al. (2016) substantially decreased in accuracy and were therefore replaced by a tailor made formula, based on manual measurements of ship hulls in Galapagos (Appendix S2). We were able to compute WSA values for 94% of all active ships in the GMR in 2016.

| Network construction
We created a network based on all available vessel tracks, consisting of geo-referenced nodes (anchorages) and edges (ship routes between these nodes). The edges were directed, that is, representing a ship route from one anchorage to another.
For each edge, we calculated two types of weights: averages of monthly aggregates of (a) ship numbers and (b) WSA. For each node, we computed unweighted in-and out-degrees (Hanneman & Riddle, 2005;Nieminen, 1974), describing the sum of all incoming and outgoing edges, respectively, using the igraph package (Csardi & Nepusz, 2006) in R. Furthermore, we calculated the weighted inand out-degrees for each node, representing the sum of all in-and outgoing flows (edge weights), respectively.

| Dispersal analysis
We considered all edges and nodes representing at least one vessel, excluding edges leading towards Outside (to avoid feedback loops), resulting in 937 edges and 106 nodes. We developed a numerical model to estimate the dispersal of potential arriving marine NIS caused by boat traffic in the GMR. A formula inspired by Seebens et al. (2013) served as a basis to calculate the risk of a species being dispersed to a node, and the capacity of that node to further spread species, which we termed dispersal signal (DS). We used the average of monthly weighted WSA edges to parameterize the amount of DS that is transferred from one node to another.
In the first iteration of the model, we assigned DS values to all receiver nodes from edges coming from Outside, where DS R is the representation of the received DS from Outside at the receiving node R, E SR all edges between the source (at iteration one the source is Outside) and the receiver, and W E the WSA at the edge E.
The constant prevents the DS from being inflated by the exponential component without changing the actual model findings (Appendix S3).
We set to 0.0001, which restricts results of the term × W E between 0 and 2. In the following iterations, all previously affected nodes are treated as source nodes that convey a certain amount of DS to the receiver nodes of all of their outgoing edges.
The extent of the conveyed DS depends on the current amount of DS at the source node DS S and the WSA value W E of the edge from the respective source to the receiver E. The new DS value for the receiver DS Rn+1 is the sum of the previous DS (DS Rn ) and the sum of conveyed DS from all incoming edges to the receiver. To simulate the dispersal risk within the GMR over several months, we excluded all edges originating During all iterations, we recorded the following parameters: (a) percentage of affected nodes and edges (visited by at least one ship), (b) accumulated DS of all nodes across the network and (c) the iteration when the network was saturated (all nodes and edges that could be reached were affected). We tested management scenarios by sequentially removing model (network) elements, that is, edges, nodes and vessel types, simulating that these elements did not contribute to the dispersal of species. This method has been proposed for the analysis of species dispersal and the spread of NIS (Keitt et al., 1997;Lookingbill et al., 2010). In total, we performed 937 edge removal, 106 node removal and nine vessel type removal scenarios.
We compared the recorded variables (a), (b) and (c) of each removal scenario to the default scenario (including all model elements) to evaluate the effectiveness of the removal.

| Most relevant nodes and edges
We selected the most relevant nodes and edges of the network based on the following criteria and will reference these as 'relevant nodes' and 'relevant edges' from here on. Relevant nodes were identified by selecting the 10 highest ranking nodes of each of the fol-

| Marine traffic in the GMR
We differentiated 352 vessels for the study period, which reflected a total WSA of 72,038 m 2 . The most abundant ship class was recreational yachts composed of small-to medium-sized private international vessels, which peaked in March and April however, they represented almost 50% of the entire WSA in the GMR in the same month.

| Node characteristics
In the constructed network, we identified 106 nodes. In-and outdegrees of nodes were similar to each other due to this network being based on ship movements, that is, in most cases the same number of ships that arrived at an anchorage also departed from it.
Out of this reason we only report results related to in-degrees.
In general, port nodes dominated the highest ranks of all indegrees and similar nodes dominated the first 7 ranks of unweighted

| Edge characteristics
The constructed network encompassed 937 edges, of which most were characterized by very few ships and WSA (Appendix S5). Of the 10 edges with the highest number of vessels passing, 6 either started or ended in a port. Only a small fraction (1%) was used by more than 15 ships per month. The majority of edges was frequented by less than 5 (92%) ships or between 5 and 15 (7%) per month ( Figure 4b).
Wetted surface area weights showed similar distribution, as 95% of edge weights were relatively low (below 5,000 m 2 ) and a large portion of the network's WSA concentrated on 1% of the edges (above 10,000 m 2 ). In contrast to weighted edges representing ship numbers, only two edges involved port nodes among the 10 highest-ranking WSA edges. However, among the list of 25 relevant edges (Appendix S5), port nodes were the most abundant, with seven edges linking to Port Baltra and Port Santa Cruz, respectively.
Besides ports, popular marine tourism sites dominated the list of the 25 relevant edges.

| Management scenarios (1): Edge removals
Three edge removal scenarios reduced the number of affected nodes in the second iteration to below 75%, representing the strongest deceleration on the increase of affected nodes. All of these three edges started in Outside and ended in a port (scenario numbers 1-3 in Figure 5a). The overall dispersal capacity for fouling species was reduced most efficiently by four edge removals, indicated by a decrease of cumulative DS down to about 80% or more (scenario numbers 2, and 4-6 in Figure 5a). In general, among the 10 edge removals that ranked highest at reducing cumulative DS, most edges either started or ended at a port (Appendix S5).   Isabela (scenario number 9 in Figure 5b). Overall DS was most efficiently reduced by removing the three ports, Port Baltra (28%), Port Santa Cruz (34%) and Port San Cristobal (44%), followed by the popular visitor site Seymour Norte (63%; Appendix S5).

| Management scenarios (3): Ship type removals
The removal of the ship type passenger had by far the strongest effect on the network, as it reduced the number of affected nodes at the second iteration to 69% and cumulative DS at the last iteration to less than 1% (scenario number 12 in Figure 5c). In comparison to passenger vessels, other ship type removals showed weak or no effects on the dispersal speed. With regard to the overall dispersal capacity, the removal of the classes tanker, general cargo and recreational vessels reduced cumulative DS to 60%, 68% and 78%, respectively (scenario numbers 13-15 in Figure 5c).

| Translating network characteristics into management implications
In general, highly connected nodes (i.e. those with a high in-degree) act as dispersal hubs, facilitating a rapid spread of NIS within the network. Therefore, management measures would be most efficient aiming to reduce the connectivity towards and between these dispersal hubs, for example, by implementing hull fouling controls along these ship routes. In the case of Galapagos, port nodes represent the main dispersal hubs of NIS, with Port Santa Cruz ranking highest in terms of unweighted and weighted in-degrees (edges, ships and WSA). All port nodes are highly interconnected, as edges between them hold the highest number of ships, and therefore are responsible for most of the network's connectivity (Figure 4b). This feature is strongest for Port Isabela, which has only few connecting edges, but very high numbers of visiting ships, almost entirely coming from other ports.
On a smaller scale, each port represents an important hub for nodes in its proximity, which are, in some cases, exclusively visited from the respective port. Examples are the edges between Port Baltra and the tourism sites Seymour Norte, Caleta Tortuga Negra and Playa Las Bachas. In addition, the dive spot Seymour Norte is so frequented that it represents a hub on its own, having as many The potential risk of introducing fouling species is a function of the WSA and number of ships, which is why these characteristics are meaningful when planning NIS prevention measures. In the GMR, a ship fuel station is located at Port Baltra, which is why many large vessels, for example, cruisers and tankers, visit this location. In contrast, most ships at Port Isabela are smaller, day-tour boats. Implementing hull controls at locations with large vessels, such as Port Baltra, would allow to cover large amounts of WSA while checking few vessels. On the other hand, locations with many smaller vessels bear a higher risk of accumulating species from various locations and also dispersing them further. The effort to control all vessel hulls would be very high, which is why performing harbour monitoring campaigns might be preferable at nodes such as Port Isabella.
Some non-port sites represented similarly large values for the number of visiting vessels and WSAs as ports due to intense traffic of passenger vessels. In fact, some edges between popular tourism sites receive higher numbers of visiting vessels than edges connected to ports (Figure 4b). Recreational yachts are equally spread across the archipelago as passenger vessels (Figure 6b), but, on average, occur in fewer numbers and are smaller than most passenger vessels (Figure 3a). Therefore, they contribute much less to the overall WSA and connectivity than passenger vessels in the GMR ( Figure 3b).

| Dispersal analysis: Disentangling the risk of NIS spreading
The dispersal model incorporates temporal marine traffic dynamics due to the underlying network being based on individual vessel movements and considers accumulated DS at source nodes. Therefore, the model accounts for varying risks at previous visited nodes quantitatively, but is missing a qualitative differentiation between the type of nodes, which might be relevant, since vessels originating from a port bear a higher risk of carrying NIS than those coming from a remote visitor site. Results of our model comply with general logic, meaning that sites with high numbers of connected edges also rank high at accumulated DS (Appendix S4). A model validation with information about NIS distributions was impossible to perform because of missing data. Standardized monitoring, for example, by settlement plates, at multiple port and non-port sites, as well as samples directly from boat hulls would be needed to validate the amounts different sites and vessel groups contribute to the overall dispersal.

| Management scenarios (1): Edge removals
The reduction of overall dispersal capacity by removing edges does not only depend on edge weights (i.e. number of passing ships and WSA values), but also on the intensity of marine traffic at the edge's arriving node, meaning the sum of all weighted and unweighted in-degrees at that node. When edges with high ship number and WSA weights were removed, the effect for reducing overall dispersal capacity was substantially lower compared to the removal of edges with lower edge weights, but therefore ending or starting in a port (Appendix S5). Ships arriving to the GMR bear the highest dispersal risk, as reflected by the strong effects of removals of edges starting in Outside on reducing the dispersal capacity in the archipelago (Figure 7). Hence, when aiming at reducing the risk of NIS spreading within the archipelago, control and enforcement of rules is most important in nodes receiving entering vessels.

| Management scenarios (2): Node removals
Node removal scenarios were more efficient at reducing the overall dispersal capacity than the removal of single edges because they reflect the summed removal of all edges connected to the respective node. Removals of ports reduced the dispersal capacity down to 28%-44% underlining their function as major dispersal hubs in the GMR (Figure 5b). The removal of Port Isabela was less efficient than the other ports because of its low values of connecting edges and WSA (Figure 7, Appendix S5). Although major port removals lead to the strongest decrease in dispersal capacity, the removals of some popular tourism sites achieve similar reductions (63%-77%; Figure 5b). This is due to cruise ships (in the passenger vessel group) travelling to the same visitor sites all year and therefore concentrating the amount of WSA to a few nodes.
These results are highly relevant for management efforts aiming to distribute tourism activity more equally across the archipelago and hence reducing impacts on frequently visited sites. We suggest establishing regular monitoring campaigns at these visitor sites to be able to respond quickly, since fast prevention measures are crucial to prevent further spreading of marine NIS (Ojaveer et al., 2015).

| Management scenarios (3): Ship classes
Removing specific ship classes represented the most unselective scenarios because they remove a whole group of vessels rather than individual network elements. Passenger vessels are the most common group of vessels in the GMR throughout the year and move across the entire archipelago (Figures 3 and 6; Appendix S6).
Thus, their removal has a tremendous effect on the dispersal capacity, reducing it to less than 1% and underlining their function as the main fouling species distributing vessel group in the GMR (Figures 5c and 7). In comparison, the removals of cargo vessels (namely tanker, general cargo and container vessels) ranked much lower at reducing dispersal capacity, since they do not visit many different nodes, but principally engage with ports only. However, due to their comparatively larger individual sizes, they contribute a substantial part of total WSA to the study area (Figure 2). This result underlines the importance of treating vessels individually in models depending on their characteristics (Floerl & Inglis, 2005;Verling et al., 2005).

| Reflecting on current biosecurity measures in Galapagos
On a global level, there is no binding regulation for hull husbandry.
However, there are some cases of regional hull management in New Zealand and Galapagos MPI, 2014). In these two examples, vessel hulls are inspected for potential biofouling upon their arrival to the respective island group. Furthermore, the Galapagos National Park Directorate (DPNG) prohibits the discharge of ballast water within the GMR and therefore eliminates this pathway. Vessels entering the GMR must navigate to one of the main ports in the archipelago, as it is not permitted to anchor anywhere but in an official port. On arrival, local authorities including the Galapagos Biosecurity Agency (ABG), DPNG, Galapagos Governing Council (CGREG), Ecuadorian Navy, Ministry of Transport (MTOP) and the Health Ministry (MSP) inspect all vessels. If a vessel is found to be transporting any kind of fouling species, it is asked to exit the GMR to be cleaned and then return for a second inspection (LOREG, 2015). Thus, the risk of introducing fouling NIS associated with cargo vessels and international recreational yachts should be substantially reduced, given the fact they only stay for a couple of days to weeks in the GMR. However, these regulations might not be able to completely prevent NIS introduction. For example, underwater hull inspections are compromised by environmental conditions, such as poor visibility and swell. In addition, hull inspections by divers are very effort-intensive, and inspections by remote operating vehicles (ROVs) are prone to overlook fouling species (Peters et al., 2019;Zabin et al., 2018). One reason for this are the so-called niche areas, which are cryptic spots on hulls, for example, rudders, propellers or sea chests, offering ideal habitats for fouling species whilst being difficult to inspect (Coutts & Taylor, 2004).
In contrast to cargo vessels and recreational yachts, domestic Galapagos vessels (namely passenger, fishing and patrol vessels) operate the whole year round in the GMR and barely leave the archipelago, which is why they bear the highest risk of spreading NIS from ports to more remote visitor sites. The largest part of the Galapagos fishing fleet is composed out of fibre reinforced plastic boats, not equipped with AIS transceivers, and therefore have not been part of our analysis. Furthermore, there are only few patrol vessels, which is why the largest part of the Galapagos fleet analysed in this study is composed of passenger vessels (Figure 3). Vessels registered in Galapagos are required to conduct a dry maintenance at least once every 2 years by the DPNG and the Ecuadorian Navy (DPNG, 2020).
In the case of larger boats, such as live-aboard cruisers, dry maintenance facilities at the islands are insufficient, causing them to visit facilities at the mainland of Ecuador. Despite these regulations, there are no obligatory hull fouling inspections for vessels operating within the archipelago other than for the boats that re-enter the GMR. Therefore, management regulations could be adjusted to mitigate the spread of marine NIS within the GMR.
Locations strongly frequented by vessels and visitors are subject to higher risks of receiving NIS but also to high amounts of disturbance, such as trash, pollution or noise. This disturbance might lower the biological resistance of communities and therefore favour the establishment of new NIS (Colautti et al., 2006). Furthermore, concentrating the propagule pressure (exposure to potential NIS) at few sites increases the chance for successful NIS establishment (Colautti et al., 2006;Seebens et al., 2019). So far, the DPNG regulates the maximum number of concurrent vessels or visitors at tourism sites, but does not apply weekly or monthly thresholds (Reck et al., 2010).
To avoid accumulations of vessels over time, we recommend tourism regulations focusing on the identified visitor sites with high numbers of passenger vessels (Appendix S4) and implement measures for the distribution of touristic activities more equally throughout the archipelago. However, this redistribution has to be managed carefully,  I G U R E 7 List of the 15 most efficient removal scenarios in terms of decreasing the overall dispersal capacity based on results of the dispersal model (compare Figure 5). The y-axis shows the removed element (edge, node or ship type) and the x-axis the relative reduced dispersal capacity in comparison to the default model (including all model elements) in percentage [Colour figure can be viewed at wileyonlinelibrary.com] et al., 2016). Although most non-indigenous fouling species in Galapagos are invertebrates  and thus difficult to identify, citizen science programmes could be steered towards a few prominent species, such as the sessile bryozoans Bugula neritina and Watersipura subtorquata. In contrast to many other fouling species, these examples are known to grow on antifouling paint and therefore are likely to be dispersed to popular marine visitor sites (Floerl & Inglis, 2005). Moreover, species watch cards and a reporting system for target species could be set up for naturalist guides and visitors. Newcomer et al. (2019) have shown that photographs taken by tourists have the potential to be used for marine invertebrate identifications, although they do not capture as many species as identifications of fouling communities in the laboratory. In Galapagos, the DiveStat program (Moity et al., 2019) has started collecting information about species sightings from tourists and park rangers, and could be complemented as to include characteristic NIS.

| CON CLUS IONS
New methodologies to analyse NIS pathways and develop risk assessments are urgently needed in an increasingly connected world.
The presented method allows to study and assess NIS dispersal risk for any setting in which ship position data are available, which is especially useful in places where ecological information about marine NIS is scarce. However, the model could be improved by weighting dispersal capacities of vessels depending on their history of previous anchorages, as well as basing the marine traffic network on actual monthly data instead of averaged monthly values. Hull controls are rarely gapless, and in Galapagos, in those cases in which marine NIS slip through the controls, there is a high risk of species spreading rapidly due to the dense shipping network. Based on our findings, we suggest to implement the following biosecurity measures in Galapagos: (a) additional hull controls at ports; (b) regular marine NIS monitoring at popular marine visitor sites and (c) adjustments of hull husbandry regulations and controls of passenger vessels, which represent the main vector for fouling NIS dispersal in Galapagos. These measures would complement existing regulations of local authorities, which focus on arriving ships to the GMR.

ACK N OWLED G EM ENTS
We thank the Galapagos National Park for the provision of the AIS data, as well as the Galapagos Biosecurity Agency (ABG) for infor-

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used to perform the described analysis contain confidential information about private vessels and their movements, which is why, according to the Galapagos National Park Directorate, it cannot be made publicly available.