Modeled larval connectivity of a multi-species reef fish and invertebrate assemblage off the coast of Moloka‘i, Hawai‘i

We use a novel individual-based model (IBM) to simulate larval dispersal around the island of Moloka‘i in the Hawaiian Archipelago. Our model uses ocean current output from the Massachusetts Institute of Technology general circulation model (MITgcm) as well as biological data on four invertebrate and seven fish species of management relevance to produce connectivity maps among sites around the island of Moloka‘i. These 11 species span the range of life history characteristics of Hawaiian coral reef species and show different spatial and temporal patterns of connectivity as a result. As expected, the longer the pelagic larval duration (PLD), the greater the proportion of larvae that disperse longer distances, but regardless of PLD (3–270 d) most successful dispersal occurs either over short distances within an island (<30 km) or to adjacent islands (50–125 km). Again, regardless of PLD, around the island of Moloka‘i, connectivity tends to be greatest among sites along the same coastline and exchange between northward, southward, eastward and westward-facing shores is limited. Using a graph-theoretic approach to visualize the data, we highlight that the eastern side of the island tends to show the greatest out-degree and betweenness centrality, which indicate important larval sources and connectivity pathways for the rest of the island. The marine protected area surrounding Kalaupapa National Historical Park emerges as a potential source for between-island larval connections, and the west coast of the Park is one of the few regions on Moloka‘i that acts as a net larval source across all species. Using this IBM and visualization approach reveals patterns of exchange between habitat regions and highlights critical larval sources and multi-generational pathways to indicate priority areas for marine resource managers.

239 generational dispersal corridors. Graph theory also provides a powerful tool for spatial 240 visualization, allowing for rapid, intuitive communication of connectivity results to researchers, 241 managers, and the public alike. This type of analysis can be used to model pairwise relationships 242 between spatial data points by breaking down individual-based output into a series of nodes 243 (habitat sites) and edges (directed connections between habitat sites). We then used these nodes 244 and edges to examine the relative importance of each site and dispersal pathway to the greater 245 pattern of connectivity around Moloka'i, as well as differences in connectivity patterns between 246 species Holstein, Paris, & Mumby, 2014). We used the R package 'igraph' to 247 examine several measures of within-island connectivity (Csardi & Nepusz, 2006). Edge density, 248 or the proportion of realized edges out of all possible edges, is a multi-site measure of 249 connectivity. Areas with a higher edge density have more direct connections between habitat 250 sites, and thus are more strongly connected. We measured edge density along and between the 251 north, south, east, and west coasts of Molokaʻi to examine possible population structure and 252 degree of exchange among the marine resources of local communities. 253 The distribution of shortest path length is also informative for comparing overall 254 connectivity. In graph theory, a shortest path is the minimum number of steps needed to connect 255 two sites. For example, two sites that exchange larvae in either direction are connected by a 256 shortest path of one, whereas if they both share larvae with an intermediate site but not with each 257 other, they are connected by a shortest path of two. In a biological context, shortest path can 258 correspond to number of generations needed for exchange: sites with a shortest path of two 259 require two generations to make a connection. Average shortest path, therefore, is a descriptive 260 statistic to estimate connectivity of a network. If two sites are unconnected, it is possible to have 261 infinite-length shortest paths; here, these infinite values were noted but not included in final 262 analyses. 263 Networks can also be broken in connected components (Csardi & Nepusz, 2006). A 264 weakly connected component (WCC) is a subgraph in which all nodes are not reachable by other 265 nodes. A network split into multiple WCCs indicates separate populations that do not exchange 266 any individuals, and a large number of WCCs indicates a low degree of island-wide connectivity. 267 A strongly connected component (SCC) is a subgraph in which all nodes are directly connected 268 and indicates a high degree of connectivity. A region with many small SCCs can indicate high 269 local connectivity but low island-wide connectivity. Furthermore, component analysis can 270 identify cut nodes, or nodes that, if removed, break a network into multiple WCCs. Pinpointing 271 these cut nodes can identify potential important sites for preserving a population's connectivity, 272 and could inform predictions about the impact of site loss (e.g. a large-scale coral bleaching 273 event) on overall connectivity. 274 On a regional scale, it is important to note which sites are exporting larvae to, or 275 importing larvae from, other sites. To this end, we examined in-degree and out-degree for each 276 region. In-degree refers to the number of inward-directed edges to a specific node, or how many 277 other sites provide larvae into site 'A'. Out-degree refers to the number of outward-directed 278 edges from a specific node, or how many sites receive larvae from site 'A'. Habitat sites with a 279 high out-degree seed a large number of other sites, and indicate potentially important larval 280 sources, while habitat sites with a low in-degree rely on a limited number of larval sources and 281 may therefore be dependent on connections with these few other sites to maintain population 282 size. Finally, betweenness centrality (BC) refers to the number of shortest paths that pass through 283 a given node, and may therefore indicate connectivity pathways or 'chokepoints' that are 284 important to overall connectivity on a multigenerational timescale. BC was weighted with the 285 proportion of dispersal as described in the preceding section. We calculated in-degree, out-286 degree, and weighted betweenness centrality for each region in the network for each species. 287 As with the source-sink index, we did not include sites on islands other than Moloka'i in 288 our calculations of edge density, shortest paths, connected components, cut nodes, in-and out-289 degree, or betweenness centrality in order to focus on within-island patterns of connectivity.

I) Effects of biological parameters on fine-scale connectivity patterns 292
The species-specific parameters that were available to parameterize the dispersal models 293 substantially influenced final output (Fig. 2). The proportion of successful settlers (either to 294 Moloka'i or to neighboring islands) varied widely by species, from 2% (Panulirus spp.) to 25% 295 (Cellana spp.). Minimum pelagic duration and settlement success were negatively correlated 296 (e.g. an estimated -0.79 Pearson correlation coefficient). Species modeled with batch spawning at 297 a specific moon phase and/or time of day (Cellana spp., P. meandrina, and C. ignoblis) 298 displayed slightly higher settlement success than similar species modeled with constant 299 spawning over specific months. On a smaller scale, we also examined average site-scale local 300 retention, comparing only retention to the spawning site versus other sites on Moloka'i ( Fig. 2). 301 Local retention was lowest for Caranx spp. (<1%) and highest for O. cyanea and P. sexfilis 302 (8.1% and 10%, respectively). 303 We measured network-wide connectivity via distribution of shortest paths, or the 304 minimum number of steps between a given two nodes in a network, only including sites on 305 Moloka'i ( Fig. 2). O. cyanea and P. sexfilis showed the smallest shortest paths overall, meaning 306 that on average, it would take fewer generations for these species to demographically bridge any 307 given pair of sites. Using maximum shortest path, it could take these species three generations at 308 most to connect sites. Cellana spp. and P. meandrina, by comparison, could take as many as five 309 generations. Other medium-and long-dispersing species showed relatively equivalent shortest-310 path distributions, with trevally species showing the highest mean path length and therefore the 311 lowest island-scale connectivity. 312 The number and size of weakly-connected and strongly-connected components in a 313 network is also an informative measure of connectivity (Fig. 2). No species in our study group 314 was broken into multiple weakly-connected components; however, there were species-specific 315 patterns of strongly connected sites. O. cyanea and P. sexfilis were the most strongly connected, 316 with all sites in the network falling into a single SCC. Cellana spp. and P. meandrina each had 317 approximately 60% of sites included in a SCC, but both show fragmentation with 7 and 6 SCCs 318 respectively, ranging in size from 2-22 sites. This SCC pattern suggests low global connectivity 319 but high local connectivity for these species. Medium and long dispersers showed larger 320 connected components; 70% of parrotfish sites fell within two SCCs; 40% of P. porphyreus sites 321 fell within two SCCs; 70% of C. strigosus sites, 55% of C. melampygus sites, and 40% of 322 Panulirus sites fell within a single SCC. In contrast, only 26% of C. ignoblis sites fell within a 323 single SCC. It is also important to note that the lower connectivity scores observed in long-324 dispersing species likely reflect a larger scale of connectivity. Species with a shorter PLD are 325 highly connected at reef and island levels but may show weaker connections between islands. 326 Species with a longer PLD, such as trevally or spiny lobster, are likely more highly connected at 327 inter-island scales which reflects the lower connectivity scores per island shown here.
374 characterized by discrete settlement pulses, whereas other species showed settlement over a 375 broader period of time. Some species also showed distinctive patterns of settlement to other 376 islands; our model suggests specific windows when long-distance dispersal is possible, as well as 377 times of year when local retention is maximized (Fig. 5).
378 III) Regional patterns of connectivity in Molokaʻi coastal waters 379 Within Moloka'i, our model predicts that coast-specific population structure is likely; 380 averaged across all species, 84% of individuals settled back to the same coast on which they 381 were spawned rather than a different coast on Molokaʻi. Excluding connections with a very low 382 proportion of larvae (<0.1% of total larvae of that species that settled to Moloka'i), we found that 383 the proportion of coast-scale local retention was generally higher than dispersal to another coast, 384 with the exception of the west coast ( Fig 6A). The north and south coasts had a high degree of 385 local retention in every species except for the long-dispersing Panulirus spp., and the east coast 386 also had high local retention overall. Between coasts, a high proportion of larvae that spawned 387 on the west coast settled on the north coast, and a lesser amount of larvae were exchanged from 388 the east to south and from the north to east. With a few species-specific exceptions, larval 389 exchange between other coasts of Moloka'i was negligible. 390 We also calculated edge density, including all connections between coasts on Moloka'i 391 regardless of settlement proportion (Fig. 6B). The eastern coast was particularly well-connected, 392 with an edge density between 0.14 and 0.44, depending on species. The southern shore showed 393 high edge density for short and medium dispersers (0.16-0.39) but low for long dispersers 394 (<0.005). The north shore also showed relatively high edge density (0.20 on average), although 395 these values were smaller for long dispersers. The west coast showed very low edge density, 396 with the exceptions of O. cyanea (0.37) and P. sexfilis (0.13). Virtually all networks that 397 included two coasts showed lower edge density. One exception was the east/south shore 398 network, which had an edge density of 0.10-0.65 except for Cellana spp. Across species, edge 399 density between the south and west coasts was 0.12 on average, and between the east and west 400 coasts was 0.04 on average. Edge density between north and south coasts was particularly low 401 for all species (<0.05), a divide that was especially distinct in Cellana spp. and P. meandrina, 402 which showed zero realized connections between these coasts. Although northern and southern 403 populations are potentially weakly connected by sites along the eastern (P. meandrina) or 404 western (Cellana spp.) shores, our model predicts very little, if any, demographic connectivity. 405 To explore patterns of connectivity on a finer scale, we pooled sites into regions (as 406 defined in Fig. 1) in order to analyze relationships between these regions. Arranging model 407 output into node-edge networks clarified pathways and regions of note, and revealed several 408 patterns which did not follow simple predictions based on PLD (Fig. 7). Cellana spp. and P. 409 meandrina showed the most fragmentation, with several SCCs and low connectivity between 410 coasts. Connectivity was highest in O. cyanea and P. sexfilis, which had a single SCC containing 411 all regions. Medium and long dispersers generally showed fewer strongly connected regions on 412 the south shore than the north shore, with the exception of C. strigosus. P. porphyreus showed 413 more strongly connected regions east of Kalaupapa but lower connectivity on the western half of 414 the island. 415 Region-level networks showed both species-specific and species-wide patterns of 416 connectivity (Fig. 8). With a few exceptions, sites along the eastern coast -notably, Cape Halawa 417 and Pauwalu Harbor -showed relatively high betweenness centrality, and may therefore act as 418 multigenerational pathways between north-shore and south-shore populations. In Cellana spp., 419 Leinapapio Point and Mokio Point had the highest BC, while in high-connectivity O. cyanea and 466 roles, even in a small MPA such as Kalaupapa. Across species, the east coast of Kalaupapa 467 showed a significantly higher betweenness centrality than the west (p=0.028), while the west 468 coast of Kalauapapa showed a significantly higher source-sink index than the east (p=2.63e-9).
469 Discussion: 470 I) Effects of biological and physical parameters on connectivity 471 We incorporated the distribution of suitable habitat, variable reproduction, variable PLD, 472 and ontogenetic changes in swimming ability and empirical vertical distributions of larvae into 473 our model to increase biological realism, and assess how such traits impact predictions of larval 474 dispersal. The Wong-Ala et al. (2018) IBM provides a highly flexible model framework that can 475 easily be modified to incorporate either additional species-specific data or entirely new 476 biological traits. In this study, we included specific spawning seasons for all species, as well as 477 spawning by moon phase for Cellana spp., P. meandrina, and C. ignoblis because such data was 478 available for these species. It proved difficult to obtain the necessary biological information to 479 parameterize the model, but as more data about life history and larval behavior become available, 480 such information can be easily added for these species and others. Some potential additions to 481 future iterations of the model might include density of reproductive-age adults within each 482 habitat patch, temperature-dependent pelagic larval duration (Houde, 1989 In this study, we have demonstrated that patterns of fine-scale connectivity around 487 Moloka'i are largely species-specific and can vary with life history traits, even in species with 488 identical pelagic larval duration. For example, the parrotfish S. rubroviolaceus and C. 489 perspicillatus show greater connectivity along the northern coast, while the goatfish P. 490 porphyreus shows higher connectivity along the eastern half of the island. These species have 491 similar PLD windows, but vary in dispersal depth and spawning season. Spawning season and 492 timing altered patterns of inter-island dispersal (Fig. 5) as well as overall settlement success, 493 which was slightly higher in species that spawned by moon phase (Fig. 2). While maximum PLD 494 did appear play a role in the probability of rare long-distance dispersal, minimum PLD appears to 495 be the main driver of average dispersal distance (Fig. 2). Overall, species with a shorter 496 minimum PLD had higher settlement success, shorter mean dispersal distance, higher local 497 retention, and higher local connectivity as measured by the amount and size of strongly 498 connected components.

499
The interaction of biological and oceanographic factors also influenced connectivity 500 patterns. Because mesoscale current patterns can vary substantially over the course of the year, 501 the timing of spawning for certain species may be critical for estimating settlement (Wren et al., 502 2016; Wong-Ala et al. 2018). Intermittent ocean processes may influence the probability of local 503 retention versus long-distance dispersal; a large proportion of larvae settled to O'ahu, which is 504 somewhat surprising given that in order to settle from Moloka'i to O'ahu, larvae must cross the 505 Kaiwi Channel (approx. 40 km). However, the intermittent presence of mesoscale gyres may act 506 as a stabilizing pathway across the channel, sweeping larvae up either the windward or leeward 507 coast of O'ahu depending on spawning site. Likewise, in our model long-distance dispersal to 508 Hawai'i Island was possible at certain times of the year due to a gyre to the north of Maui; larvae 509 were transported from Kalaupapa to this gyre, where they were carried to the northeast shore of 510 Hawai'i (Fig. S6). Preliminary analysis also suggests that distribution of larval depth influenced 511 edge directionality and size of connected components (Fig. 7); surface currents are variable and 512 primarily wind-driven, giving positively-buoyant larvae different patterns of dispersal than 513 species that disperse deeper in the water column (Fig. S7).

II) Model limitations and future perspectives 515
Our findings have several caveats. Because fine-scale density estimates are not available 516 for our species of interest around Moloka'i, we assumed that fecundity is equivalent at all sites. 517 This simplification may lead us to under-or over-estimate the strength of connections between 518 sites. Lack of adequate data also necessitated estimation or extrapolation from congener 519 information for larval traits such as larval dispersal depth and PLD. Since it is difficult if not 520 impossible to identify larvae to the species level without genetic analysis, we used genus-level 521 larval distribution data (Boehlert & Mundy 1996), or lacking that, an estimate of 50-100m as a 522 depth layer that is generally more enriched with larvae (Boehlert et al., 1992, Wren & 523 Kobayashi, 2016). We also estimated PLD in several cases using congener-level data (see Table  524 1). While specificity is ideal for making informed management decisions about a certain species, 525 past sensitivity analysis has shown that variation in PLD length does not greatly impact patterns 526 of dispersal in species with a PLD of >40 days . 527 Although our MITgcm current model shows annual consistency, it only spans two and a 528 half years chosen as neutral state 'average' ocean conditions. It does not span any El Niño or La 529 Niña (ENSO) events, which cause wide-scale sea-surface temperature anomalies and may 530 therefore affect patterns of connectivity during these years. El Niño can have a particularly 531 strong impact on coral reproduction, since the warm currents associated with these events can 532 lead to severe temperature stress (Glynn & D'Croz, 1990;Wood et al., 2016). While there has 533 been little study to date on the effects of ENSO on fine-scale connectivity, previous work has 534 demonstrated increased variability during these events. For example, Wood et al. (2016) showed 535 a decrease in eastward Pacific dispersal during El Niño years, but an increase in westward 536 dispersal, and Treml et al. (2008) showed unique connections in the West Pacific as well as an 537 increase in connectivity during El Niño. While these effects are difficult to predict, especially at 538 such a small scale, additional model years would increase confidence in long-term connectivity 539 estimations. Additionally, with a temporal resolution of 24 hours, we could not adequately 540 address the role of tides on dispersal, and therefore did not include them in the MITgcm. 541 Storlazzi et al. (2017) showed that tidal forces did affect larval dispersal in Maui Nui, 542 underlining the importance of including both fine-scale, short-duration models and coarser-scale, 543 long-duration models in final management decisions. 544 We also limit our model's scope geographically. Our goal was to determine whether we 545 could resolve predictive patterns at this scale relevant to management. Interpretation of 546 connectivity output can be biased by spatial resolution of the ocean model, since complex coastal 547 processes can be smoothed and therefore impact larval trajectories. To limit this bias, we focused 548 mainly on coastal and regional connectivity on scales greater than the current resolution. We also 549 used the finest-scale current products available for our study area, and our results show general 550 agreement with similar studies of the region that use a coarser resolution  and 551 a finer resolution (Storlazzi et al., 2018). Also, while knowledge of island-scale connectivity is 552 important for local management, it does disregard potential connections from other islands. In 553 our calculations of edge density, betweenness centrality, we included only settlement to 554 Moloka'i, discarding exogenous sinks that would bias our analysis. Likewise, we cannot predict 555 the proportion of larvae settling to other islands that originated from Moloka'i, or the proportion 556 of larvae on Moloka'i that originated from other islands. 557 It is also important to note scale in relation to measures of connectivity; we expect that 558 long-dispersing species such as Caranx spp. and Panulirus spp. will show much higher measures 559 of connectivity when measured across the whole archipelago as opposed to a single island. The 560 cut-nodes observed in these species may not actually break up populations on a large scale due to 561 this inter-island connectivity. Nevertheless, cut-nodes in species with short-and medium-length 562 PLD may indeed mark important habitat locations, especially in terms of providing links 563 between two otherwise disconnected coasts. It may be that for certain species or certain regions, 564 stock replenishment relies on larval import from other islands, underscoring the importance of 565 MPA selection for population maintenance in the archipelago as a whole.

III) Implications for management 567
Clearly, there is no single management approach that encompasses the breadth of life 568 history and behavior differences that impact patterns of larval dispersal and connectivity (Toonen 569 et al., 2011, Holstein, Paris, & Mumby, 2014. The spatial, temporal, and species-specific 570 variability suggested by our model stresses the need for multi-scale management, specifically 571 tailored to local and regional connectivity patterns and the suite of target species. Even on such a 572 small scale, different regions around the island of Moloka'i can play very different roles in the 573 greater pattern of connectivity (Fig. 8); sites along the west coast, for example, showed fewer 574 ingoing and outgoing connections than sites on the north coast, and therefore may be more at risk 575 of isolation. Seasonal variation should also be taken into account, as mesoscale current patterns 576 (and resulting connectivity patterns) vary over the course of a year. Our model suggests species-577 specific temporal patterns of settlement (Fig. 5); even in the year-round spawner O. cyanea, local 578 retention to Moloka'i as well as settlement to O'ahu was maximized in spring and early summer, 579 while settlement to other islands mostly occurred in late summer and fall. 580 Regions that show similar network dynamics may benefit from similar management 581 strategies. Areas that act as larval sources either by proportion of larvae (high source-sink index) 582 or number of sites (high out-degree) should receive management consideration. On Moloka'i, 583 across all species in our study, these sources fell mostly on the northern and eastern coasts. 584 Maintenance of these areas is especially important for downstream areas that depend on 585 upstream populations for a source of larvae, such as those with a low source-sink index, low in-586 degree, and/or low local retention. Across species, regions with the highest betweenness 587 centrality scores fell mainly in the northeast (Cape Halawa and Pauwalu Harbor). These areas 588 should receive consideration as potentially important intergenerational pathways, particularly as 589 a means of connecting north-coast and south-coast populations, which showed a lack of 590 connectivity both in total number of connections (edge density) and proportion of larvae. Both of 591 these connectivity measures were included because edge density includes all connections, even 592 those with a very small proportion of larvae, and may therefore include rare dispersal events that 593 are of little relevance to managers. Additionally, edge density comparisons between networks 594 should be viewed with the caveat that these networks do not necessarily have the same number 595 of nodes. Nevertheless, both edge density and proportion show very similar patterns, and include 596 both demographically-relevant common connections as well as rare connections that could 597 influence genetic connectivity. 598 Management that seeks to establish a resilient network of spatially managed areas should 599 also consider the preservation of both weakly-connected and strongly-connected components, as 600 removal of key cut-nodes (Fig. S5) breaks up a network. Sites within a SCC have more direct Target taxa selected for the study, based on cultural, ecological, and/or economic importance.   Manuscript to be reviewed   Dispersal distance density kernels.
Dispersal distance is combined across species by minimum pelagic larval duration (PLD) length in days (short, medium, or long). Most short dispersers settle close to home, while few long dispersers are retained at or near their spawning sites.  Region-level summary statistics across all species.
Betweenness centrality is a measure of the number of paths that pass through a certain region; a high score suggests potentially important multi-generation connectivity pathways.
In-degree and out-degree refer to the amount of a node's incoming and outgoing connections. Betweenness centrality, in-degree, and out-degree have all been normalized to values between 0 to 1 per species. Local retention is measured as the proportion of larvae that settled back to their spawn site out of all larvae spawned at that site. Source-sink index is a measure of net export or import; negative values (blue) indicates a net larval sink, while positive values (red) indicates a net larval source. White indicates that a site is neither a strong source nor sink. Gray values for Cellana spp. denote a lack of suitable habitat sites in that particular region.