To test how biotic resistance and abiotic filters interact with propagule pressure to predict invasions and to explore the architecture of shipping networks, we selected three port regions to serve as models for distinct biogeographic regions and shipping profiles. This approach allowed us to assess invasion risk in distinct ecological and invasion contexts. We used three main criteria to select our study regions. Port regions needed to: 1) be intensively used and relevant economically to the overall region, including a high volume of international connections; 2) be located at different biogeographic regions within the marine biome to represent different ecological contexts; and 3) have distinct shipping regimes with variable propagule pressure. We selected three port regions along the northeast Pacific coastline, from subarctic latitudes to the tropics. Port regions combined all nearby commercial ports (less than 100 km distance) that would be relevant to that bioregion: 1) Panama Bay anchorage region (8° 52’ 12’’ N, 79° 30’ 0’’ W), including all arrivals and anchored ships about to cross the Panama Canal; 2) California’s San Francisco Bay port region, including all arrivals to the ports and commercial harbors at Alameda (37° 45’ 36’’ N, 122° 18’ 36’’ W), Benicia (38° 1’ 48’’ N, 122° 19’ 12’’ W), Oakland (37° 47’ 60’’ N, 122° 16’ 48’’ W), San Francisco (37° 47’ 60’’ N, -122° 25’ 12’’ W), and Redwood City (37° 30’ 0’’ N, 122° 13’ 12’’ W); and 3) Southeast Alaska port region, represented by the Port of Ketchikan (55° 20’ 33’’ N, 131° 39’ 22’’ W). These three port regions (hereafter referred as Panama, California, and Alaska respectively) were also identified consistently in the literature as potential routes of introduction for marine species (see databases curated by Fofonoff et al., 2018), as well as biogeographic regions where biotic interactions by predation have shown consistent differences (Freestone et al., 2021; López and Freestone, 2021).

Panama – Panama served as a tropical latitude port region with extremely high vessel traffic that was consistent year-round, which ensured potential high propagule pressure all year. The Panama Canal is a transoceanic aquatic corridor connecting the Pacific and Atlantic across the Isthmus of Panama. It is an important biogeographical connector and potential invasion hotspot (Ruiz et al., 2009; Muirhead et al., 2015). As expansion projects are implemented over time, the growing traffic capacity of larger vessels will likely increase invasion risk in the region (Ruiz et al., 2009, 2015; Castellanos-Galindo et al., 2020).

California – Recognized as a highly invaded estuary (Cohen and Carlton, 1998; Ruiz et al., 2011), the San Francisco Bay ports are located at a temperate latitude and are major ports of entry for goods to the United States. Continuous arrivals year-round and high arrival numbers contribute to high potential propagule pressure to this region.

Alaska – Ketchikan is a major port for passenger vessels in the Northeast Pacific facilitating tourism in southeast Alaska (McGee et al., 2006; Miller and Ruiz, 2014). In the region where it is located, Ketchikan has a higher number of arriving vessels relative to neighboring ports, with arrivals that peak in frequency during the northern hemisphere summer. The importance of this port region is growing given the emerging shipping routes in the Arctic Ocean that could potentially increase risk of invasions in this region (Verling et al., 2005).

We obtained the historical port transit profile and shipping traffic information between 2006 and 2010 for the focal port regions. Data were obtained from the National Ballast Information Clearinghouse (National Ballast Information Clearinghouse, 2021) for those ports located on coastal United States of America and provided by the Panamanian Canal authorities (ACP) for Panama. Overall number of arrivals, monthly traffic, port connections (i.e., potential donor ports) and vessel type were used to estimate contributions of shipping intensity to propagule pressure (Miller et al., 2018). While our data from Panama only includes ship transits through the Panama Canal and do not include arrivals at Pacific ports that do not also transit the Canal, Pacific anchorages represent most of the annual commercial vessel traffic in the region. Therefore, our estimate of shipping traffic to this region is likely underestimated but captures the majority of regional vessel traffic. In addition, vessel traffic data from Panama initially included all loading port stops within a route for cargo vessels. In comparison, NBIC data only includes the last port-of-call for each arrival. To ensure direct comparison to the other port regions, we estimated last port-of-call for Panama arrivals by assuming that the geographically closest port was the last stop before the final destination in Panama. We then selected the smallest geodesic distance from each stop and the Panama port region within a ship route r using package Geosphere 1.5–10 (Hijmans, 2019). Invasion risk models (P[Inv]) were computed for the Panama dataset with estimated last port-of-call and all stops, but no substantial differences were observed in the results of the model. Therefore, we focus our results on the estimated last port-of-call, which allows a more direct comparison to the other studied port regions.

We further focused on hull fouling (i.e., biofouling) as an important mechanism for accidental introduction of propagules overseas. To calculate the number of individuals in each release event (i.e., ship arrival), we estimated the number of potential biofouling organisms that could colonize the wetted surface area of commercial vessel hulls (Moser et al., 2016) and associated structures (“niche areas”; rudders, propellers, shafts etc., see Moser et al., 2017; Miller et al., 2018) in each arrival. Wetted surface area (WSA) refers to the estimated area of a ship’s hull that is submerged when loaded with maximum cargo and has important applications in maritime engineering. To estimate the potential colonization space for macrofouling organisms on the analyzed vessels, we primarily employed mean total WSA with incorporated fractions of total niche area for each vessel type (Miller et al., 2018). When total niche area data was unavailable for specific vessel types (e.g., fishing vessel), we utilized total wetted surface area as a substitute (Moser et al., 2016), but those vessel arrivals represented less than 1% of our data. Both approaches serve as proxies for propagule supply, allowing us to assess the potential threat of accidental biofouling organism transfer by ships.

The likelihood of propagules being released from a ship’s hull in a recipient region j is estimated for each ship arrival using a modified equation from (Wonham et al., 2013).

Probability of propagule arrival and release, or P_r(Supply), reflects the likelihood of a successful propagule release from a ship arrival on route r. P_r(Supply) considers a fixed baseline probability (p) of 0.00002 (as estimated by Seebens et al., 2013), that represents a single propagule failing to be released and N_c , which acts as a proxy for propagule pressure based on the colonization potential of the ship’s hull. This model therefore leverages vector-based dynamics to estimate invasion risk due to the abundance of propagules per ship arrival (Wonham et al., 2013). The core principle is that the probability of colonization increases with the successful release of propagules, which P_r(Supply) captures through the interplay of p and N_c . To estimate N_c , we calculated the number of potential macrofouling organisms that could colonize each ship’s hull (WSA) using the average area occupied by a common fouling barnacle (6.45 cm²), following established methods (Moser et al., 2016). Hull fouling treatment is not as regulated as ballast water transfer overseas (Zabin et al., 2018), so this model assumes an upper limit of propagule pressure by considering every entry event to have an equal probability of propagule arrival and release, scaled to total WSA of the various ships. This approach disregards any husbandry treatment that could reduce chances of propagule release, such as specific anti-fouling coatings or hull cleaning practices. While complete hull fouling is uncommon in modern ships (Chan et al., 2022) this model serves to capture the potential magnitude of commercial vessels to deliver propagules via biofouling across recipient ports. While hull fouling is an important vector of marine invasions, ships could also transport native species within their range. Since there is no known probability that a species that is transported via hull fouling is non-native to a recipient port, this model assumes that all introduced propagules to a port have the potential to be non-native, allowing the calculation of relative invasion risk based on overall propagule pressure.

To explore the potential for environmental conditions at the arrival port to mediate invasion risk, we obtained environmental data for all connected ports in our database from published data on global ports (Keller et al., 2011). Water salinity and minimum, maximum, and average annual water temperature were obtained for each port. We consulted open-source tools such as MarineTraffic (marinetraffic.com) and VesselFinder (https://www.vesselfinder.com/) to resolve inconsistencies in port names. We further estimated missing environmental values using publicly available data (United States Geological Survey, through waterdata.usgs.gov; Bio-ORACLE, through Tyberghein et al., 2012; Assis et al., 2018). Ports that could not be identified were excluded from the final dataset but were less than 1% of all entries combined (< 400).

To explore the contribution of environmental distance to invasion risk within the shipping network, we built weighted network maps using R package visNetwork 2.0.9 (Almende et al., 2019) and igraph (Csardi and Nepusz, 2006). This approach allowed us to visualize the over 40,000 arrivals in our combined datasets and evaluate patterns of evenness (Pielou’s evenness index J′) and diversity (Shannon diversity H′) of donor ports for each study region. Evenness and diversity metrics were computed using the vegan 2.5–6 package (Oksanen et al., 2022). Environmental distances were calculated by estimating the Euclidean distance between data matrices of donor and recipient port environmental conditions (Keller et al., 2011).

The likelihood of an organism surviving abiotic filters increases with environmental similarity between donor and recipient regions, which was modelled as a modified Gaussian function (Seebens et al., 2013).

The probability of environmental survival results from the differences in mean annual water temperature T and mean annual salinity S between donor port i and recipient region j, standardized by the width of the ecological niche σT and σS, i.e., the variance of the distributions ( Supplementary Table 1 ). The base probability α was modelled with the assumption that organisms would survive perfectly in matched environmental conditions; therefore, we set values of α equal to 1. In this model either salinity or temperature can affect mismatch to the same extent (correlation tests between salinity and temperature revealed no significant relationship, ρ = 0.2). We used mean salinity and temperature conditions to estimate environmental matching given their utility in prior modeling efforts (e.g., Seebens et al., 2013, 2016, 2018), while noting that other attributes may also contribute to environmental resistance (e.g., temperature variation, upwelling regimes, etc.). The constant e is a fundamental mathematical constant, approximately equal to 2.71828, that is widely used in exponential and logarithmic functions.

To quantify predation pressure in recipient regions, and therefore the potential for biotic resistance, we conducted a short-term predator exposure experiment on biofouling communities. Experiments were deployed at three recreational marinas in each study region (see Freestone et al., 2021 for detailed experiments and results). Median distance between study sites and focal port regions was 5.8 km. Biofouling communities represent a diverse multiphyletic assemblage of marine invertebrate taxa, including barnacles, ascidians, and bryozoans capable of colonizing artificial substrates and ship’s hulls (Godwin, 2003; Sylvester et al., 2011), and these communities also harbor many of the known coastal invertebrate invasions in the United States (Ruiz et al., 2000). To control for potential confounding factors of habitat type (i.e., substrate) and area, we used polyvinyl chloride panels (PVC, 14 x 14 x 0.95 cm) as artificial model habitat for the fouling communities. While these panels may not fully replicate the complexities of real ship hulls, they offer a well-established and comparable method for assessing fouling communities on artificial substrates (Marraffini et al., 2017). Biofouling communities developed on panels for a period of three months under reduced predation (inside a marine plastic cage, mesh size 6.35 mm) and were then exposed to ambient predation for three days (n = 5 panels/site). Biomass of the community (i.e., wet weight, g) was measured before and after the experiment (see Freestone et al., 2021 for results and Freestone et al., 2022 for publicly available data). Experiments were deployed from floating docks at 1m below the water surface, allowing us to assess the impacts of predatory fish as an important predator guild affecting biofouling on visiting vessels.

We then used the experimentally-derived measures of predation impact on fouling prey communities to obtain an estimated predation survival probability as a metric for biotic resistance for each focal port region.

Probability of predation survival P_rj(Pred) at recipient port j from a ship arrival from route r is a probability density function, similar to Equation 3, and results from biofouling community biomass differences Δ W b a before (b) and after (a) exposure to predation, standardized by the variance in weight data (σW) ( Supplementary Table 1 ). The baseline probability of an organism surviving predation is indicated by α, which we set at 1. P(Pred) uses loss of biomass from predation as an estimate of surviving predation impact and overcoming biotic resistance potential. To capture the inherent variability in predation pressure, we modeled the probability of predation (P[Pred]) using a Gaussian function. This function best approximates the distribution of community biomass differences ( Δ W b a ) in each region, but does not account for any selectivity of predators on native or non-native prey, or use of habitat refuges. Therefore, the probability of predation survival increases with lower impacts of predators on biomass of fouling communities, therefore decreasing biotic resistance.

Probability of biotic resistance was then estimated using Equation 4 for each experimental sampling unit (i.e., panel) using the observed pre- and post-exposure biomass differences, rendering 15 values for each of the three study regions. To incorporate observed variability in predation pressure within and among regions into the invasion risk model, one value of P(Pred) was randomly pulled from bootstrapped values (reps = 1000) in each iteration of the model (i.e., each ship arrival) when calculating P(Inv). Therefore, for each route r a random value was selected from the calculated P(Pred) values for that recipient region j.

This model captures the magnitude of commercial vessels contributing to the potential transfer of marine organisms and its interplay with factors that could mediate the risk of invasion, including environmental mismatch and biotic resistance from predation. Given that all three probabilities contribute to risk of invasion in a recipient port and are considered independent events, their product is what determines the relative likelihood of new invasions in the port region of interest for each individual ship arrival (Equation 1). We consider, however, that risk is cumulative through time. We therefore provide estimates of invasion risk per arrival as well as a cumulative (i.e., summed) risk for the full four-year duration of the dataset (2006–2010).

Shipping intensity varied by an order of magnitude among port regions with most arrivals in Panama during 2006–2010 ( Figure 1 ). From the overall 47,456 arrivals in the dataset, 30,218 were to Panama in comparison to 15,336 and 1,902 to California and Alaska, respectively. Most vessels arriving to Alaska were passenger ships (n = 1,786), while arrivals to both Panama and California were primarily container ships (n_Panama = 10,619 and n_California = 10,431) ( Figure 1 ). There was no consistent inter-annual variation in any of the focal port regions, but Alaska had intra-annual variation with peak arrivals occurring between May and August, coinciding with expected patterns of tourism. Panama also displayed the largest diversity of vessel types, highest total WSA (Panama: 1.644e+08 vs. California: 0.994e+08, Alaska: 4.645e+06) and had an estimated total propagule size (i.e., sum of estimated fouling barnacles given WSA; Panama: 2.549e+09) that was almost 1.7 times larger than even California (1.543e+09; Alaska: 7.201e+07).

Lower environmental distance per connection increased average risk of invasions. Despite California having a number of arrivals that was an order of magnitude higher than Alaska, in both regions the number of connections per donor port was skewed with some ports being more highly connected than others, rendering very similar evenness values (California J′= 0.520; Alaska J′ = 0.526). Furthermore, Alaska had the lowest diversity of port connections ( Figure 2 ; H′ = 1.98), and the average environmental distance of connections was two times lower than the other regions, thereby increasing the risk of invasion from environmental matching given the high number of connections to environmentally close ports ( Figure 3 ). In contrast, California had a higher diversity of connections (California H′ = 2.29) to more environmentally dissimilar regions, which moderated invasion risk. While there were shared connections among all ports, Panama displayed the highest number of unique connections and therefore the highest port diversity ( Figure 2 ; H′ = 4.34), approximately twice that of California or Alaska. Panama also had a more even distribution of connections among donor ports relative to the other regions (J′ = 0.68) in addition to greater variability of environmental distances between connections. Overall, results suggest that environmental filtering might systematically reduce invasion risk in all regions but act more strongly in California.

When considering both propagule supply (P[Supply]) and resistance (P[Env] and P[Pred]) together, our models estimated a very low mean probability of invasion per arrival across all port regions ( Figure 4 ), and this low probability stemmed from a combination of environmental and biotic resistance. The relative importance of each of these factors to invasion risk, however, differed significantly among regions. Mean probability of invasion due to propagule supply alone was at least 30% higher in Panama and California than Alaska, but variability (i.e., standard deviation) was greater in Panama. Despite high propagule pressure in some regions, environmental mismatch contributed to a substantial number of potential failures in all regions and markedly reduced invasion probability (P[Inv]) in California. Further, stronger predation in tropical Panama lowered the average probability of invasion due to biotic resistance by up to 80%. In contrast, biotic resistance had a negligible effect on invasion risk in California and Alaska. Therefore, the low overall risk of invasion per arrival across all regions can be attributed to the interplay of strong biotic resistance in Panama, low propagule supply in Alaska and consistently low environmental matching in all regions.

While average invasion risk per arrival was low in all three port regions, cumulative risk due to the sheer volume of arrivals in some ports showed markedly different patterns ( Figure 5 ). In Alaska, cumulative risk increased very steeply but reached a maximum value that was forty times lower than California and Panama, due to fewer arrivals in this region ( Figure 5A , inset). Biotic resistance reduced cumulative risk of invasions in Panama, but not in the other regions ( Figure 5 ). Interestingly, cumulative risk of invasion in California was substantially reduced by environmental matching but was as high as Panama, which had double the number of arrivals but with much greater potential for biotic resistance.