Higher-order patterns of aquatic species spread through the global shipping network

The introduction and establishment of nonindigenous species (NIS) through global ship movements poses a significant threat to marine ecosystems and economies. While ballast-vectored invasions have been partly addressed by some national policies and an international agreement regulating the concentrations of organisms in ballast water, biofouling-vectored invasions remain largely unaddressed. Development of additional efficient and cost-effective ship-borne NIS policies requires an accurate estimation of NIS spread risk from both ballast water and biofouling. We demonstrate that the first-order Markovian assumption limits accurate modeling of NIS spread risks through the global shipping network. In contrast, we show that higher-order patterns provide more accurate NIS spread risk estimates by revealing indirect pathways of NIS transfer using Species Flow Higher-Order Networks (SF-HON). Using the largest available datasets of non-indigenous species for Europe and the United States, we then compare SF-HON model predictions against those from networks that consider only first-order connections and those that consider all possible indirect connections without consideration of their significance. We show that not only SF-HONs yield more accurate NIS spread risk predictions, but there are important differences in NIS spread via the ballast and biofouling vectors. Our work provides information that policymakers can use to develop more efficient and targeted prevention strategies for ship-borne NIS spread management, especially as management of biofouling is of increasing concern.

interactions (e.g. competition, predation) come into play. Global studies of biofouling richness have found 23 that the tropics have higher richness at the regional scale [13], although higher predation rate in the 24 tropics decreases local richness more compared to subtropical or temperate regions [14]. On average, bare 25 space colonization rate (the number of new species or individuals per unit time) increases with water 26 temperature [15] and fluctuates seasonally in temperate climates but is more stable in tropical climates. 27 To predict biofouling accumulation on ships the data available to us (e.g., ship duration of stay in 28 ports, ship type, port locations), we used biofouling richness data from studies that monitored coastal 29 sites over time at temperate and tropical latitudes [15,16,13,14]. In order to compare the data from 30 these studies, which used different methods and recorded different taxonomic groups, we calculated the 31 proportion of maximum species richness at each site for each study. We then fit 3rd-order polynomial 32 regression models to tropical and temperate proportional richness data separately, as temperate species 33 richness was overall lower and predicted to show different seasonal dynamics. Biofouling richness in the 34 tropics was found to peak at 4 months (∼120 days) just above 80% maximum observed richness and then 35 declined slowly to just under 50% maximum richness over the rest of the year, until a slight uptick in the 36 last month Figure S1. The temperate richness trend line increased more slowly to a peak of just under 37 50% at 7 months (210 days) followed by a slow decline to just above 35% at the end of a year. AIC model 38 selection considering 1st-order, 2nd-order, and 3rd-order polynomial models for tropical and temperate 39 data ranked the 3rd-order polynomial model best for tropical data.  Considering both the time spent in port and the likelihood of antifouling system presence, we estimate 53 biofouling species richness on a ship at source port i as a proportion of maximum richness as follows: Where d i is the time (in days) that ship v stays at source port i, and A (t) refers to the proportion of 55 ships of a given type without an operational antifouling system. β T r1 = 1.29 × 10 −7 , β T r2 = 8.316 × 10 −5 , 56 β T r3 = 0.0149 are the tropical coefficients and β T p1 = 1.
are the temperate coefficients. Where v ij is the average voyage velocity of ship v from source port i to destination port j measured in 59 kilometers per day, and γ = 0.008. This model (displayed in Figure S2) yields a survival rate of 0.007 at 60 14 knots (or 622.2 kilometers per day), which is similar to the survival rate experimentally found in [18].

61
Total Voyage Introduction Risk We estimate the relative biofouling risk of each ship voyage p ij in the Lloyd's dataset by multiplying equations 1 (or 2) and 3. If | source port latitude | < 35.0 decimal degrees: If | source port latitude| ≥ 35.0 decimal degrees:

Caveats
As mentioned at the beginning of this section, biofouling is a complex process that is 62 dependent on many environmental, ship, and voyage factors. For the purposed of our research, we focused 63 on modeling biofouling introduction risk as a function of parameters that 1) could be predicted within 64 our shipping and environmental dataset in a global scale and 2) could be estimated with empirical data.

65
Other parameters that are known to influence biofouling introduction risk which could not be included in       Table S2: p-values obtained from the t-tests performed between SF-HON models and the two other models.
dataset in which ballast SF-HON has a lower error across all states. We further illustrate the prediction 113 of each model versus the actual introductions in Figure S7. Overall, we observe an under-estimation of 114 the risk by SF-HON predictions and over-estimation of the risk by All-Paths model predictions. It is 115 important to note that, some states/countries had fewer records resulting in higher prediction error and 116 higher variance, which is a limitation of our work.