Molecular Ecology (2014) 23, 5036–5047

doi: 10.1111/mec.12924

Ocean circulation model predicts high genetic structure observed in a long-lived pelagic developer J . M . S U N D A Y , * † I . P O P O V I C , † W . J . P A L E N , † M . G . G . F O R E M A N ‡ and M . W . H A R T † *Biodiversity Research Centre, University of British Columbia, 2212 Main Mall, Vancouver, British Columbia, Canada, †Department of Biological Sciences, Simon Fraser University, 8888 University Drive, Burnaby, British Columbia, Canada, ‡Institute of Ocean Sciences, Fisheries and Oceans Canada, 9860 West Saanich Road, Sidney, British Columbia, Canada

Abstract Understanding the movement of genes and individuals across marine seascapes is a long-standing challenge in marine ecology and can inform our understanding of local adaptation, the persistence and movement of populations, and the spatial scale of effective management. Patterns of gene flow in the ocean are often inferred based on population genetic analyses coupled with knowledge of species’ dispersive life histories. However, genetic structure is the result of time-integrated processes and may not capture present-day connectivity between populations. Here, we use a high-resolution oceanographic circulation model to predict larval dispersal along the complex coastline of western Canada that includes the transition between two well-studied zoogeographic provinces. We simulate dispersal in a benthic sea star with a 6–10 week pelagic larval phase and test predictions of this model against previously observed genetic structure including a strong phylogeographic break within the zoogeographical transition zone. We also test predictions with new genetic sampling in a site within the phylogeographic break. We find that the coupled genetic and circulation model predicts the high degree of genetic structure observed in this species, despite its long pelagic duration. High genetic structure on this complex coastline can thus be explained through ocean circulation patterns, which tend to retain passive larvae within 20–50 km of their parents, suggesting a necessity for close-knit design of Marine Protected Area networks. Keywords: larval dispersal, marine connectivity, oceanographic circulation model, population genetics Received 1 June 2013; revision received 9 September 2014; accepted 12 September 2014

Introduction Understanding patterns of gene flow across geographic landscapes represents an important intersection between the goals of ecology and evolution. The movement of genes and individuals across landscapes can provide information on the possibility of local adaptation, the persistence of populations, and the spatial scale at which management efforts will be most effective. The use of molecular markers to infer gene flow from the spatial distribution of alleles has become widespread, but is fraught with potential errors due to Correspondence: J. M. Sunday, Fax: (604) 827 5350; E-mail: [email protected]

mismatches between the evolutionary processes that may lead to genetic differentiation (e.g. genetic drift and historical vicariance) and the highly simplistic assumptions often made when interpreting differentiation as a consequence of contemporary gene flow (Whitlock & McCauley 1999; Marko & Hart 2011). For example, Wright’s Island model is widely used to calculate and interpret genetic structure using F-statistics, but this interpretation depends on the assumption of drift–migration equilibrium and attributes allele sharing between populations only to gene flow rather than to alleles shared by descent from an ancestral population (Wright 1969; Slatkin 1987). Coalescent approaches, such as the isolation-with-migration model, provide a nonequilibrium-based framework in which it is possible © 2014 John Wiley & Sons Ltd

O C E A N C I R C U L A T I O N M O D E L P R E D I C T S G E N E T I C S T R U C T U R E 5037 to distinguish the effects of common ancestry from gene flow (Nielsen & Wakeley 2001; Hey & Nielsen 2004). However, gene flow estimates from these methods are not very sensitive to recent demographic changes and tend to require data from many loci (Marko & Hart 2011). A complementary approach for using genetic structure to understand connectivity between populations is to have a priori predictions of the scale and pattern of population genetic structure based on independent ecological data (Irwin 2002). Independent information on migratory behaviour, dispersal probabilities or evidence of localized adaptation can help reinforce or discount gene flow as the cause of population genetic patterns (Irwin 2002). In marine systems, population connectivity is thought to depend on oceanographic currents combined with characteristics of dispersive life-history stages of taxa (Palumbi 2004; Levin 2006). While dispersal patterns are notoriously intractable, high-resolution oceanographic models are becoming increasingly available to provide predictions of dispersal and thus connectivity between geographic locations as a source of information to complement genetic studies (Siegel et al. 2003). Several recent studies have used such oceanographic models for direct comparison with empirical genetic data. These vary from simple qualitative comparisons of dispersal patterns with observed genetic structure (e.g. Gilg & Hilbish 2003; Baums et al. 2006; Kenchington et al. 2006), to analytical tests comparing estimated dispersal probabilities with observed pairwise genetic differentiation (Galindo et al. 2006, 2010; Dupont et al. 2007; White et al. 2010; Foster et al. 2012), or multilocus assignment tests (Schunter et al. 2011). Results of these studies have either shown good agreement with genetic studies (Galindo et al. 2006, 2010; White et al. 2010; Foster et al. 2012) or have been useful for identifying alternative or additional causes of population structure, such as the role of history, natural selection and larval behaviour (e.g. Galindo et al. 2010). Although there have been several studies quantitatively comparing oceanographic patterns to genetic structure, little work has been done in regions strongly influenced by cyclical Pleistocene glaciations where genetic structure is most readily attributed to historical vicariance or population expansions and recolonization. The role of contemporary gene flow acting on short ecological timescales in such history-dominated regions is often overlooked. Yet, in such glacially influenced coastlines, deep fjords and complex sea-floor topography can cause oceanographic features that may entrain dispersive propagules. Oceanographic structure in such regions thus offers an important null expectation against which the influence of historical events can be © 2014 John Wiley & Sons Ltd

identified from genetic patterns. Moreover, it offers critical information at the timescale needed for management applications such as the design of marine protected areas. The northwest coast of North America is a region of complex glacial history and oceanography and has attracted many population genetics studies. Cyclical glacial cover over much of the coastal region during the Pleistocene likely affected the biogeographic distributions of coastal marine species (Pielou 1991). This glacier-pocketed coastline includes large fjords, canyons and islands, which promote the formation of multiple eddies and small countercurrents (Foreman et al. 2008). In addition, two major oceanic gyres split along the coastline within this region, forming the Alaska current moving northwards, and the California current moving southwards, along with their associated coastal countercurrents (Dodimead et al. 1963). This split roughly coincides with the transition between the Oregonian and Aleutian zoogeographical provinces (Briggs 1974). Numerous comparative studies of genetic structure in this region have tended to reveal either relatively low diversity and little genetic structure throughout, or relatively high diversity with significant genetic structure among populations from Oregon to Alaska (Hellberg 1996; Arndt & Smith 1998; Rocha-Olivares & Vetter 1999; Kyle & Boulding 2000; Edmands 2001; Marko 2004; Hickerson & Cunningham 2005; Harley et al. 2006; Wilson 2006; Kelly & Eernisse 2007; Kelly & Palumbi 2010). Most strikingly, these patterns do not seem to be related to pelagic larval duration (Marko 2004; Kelly & Palumbi 2010), and several species with long pelagic durations have been found to have significant population structure across this region (Rocha-Olivares & Vetter 1999; Kyle & Boulding 2000; Keever et al. 2009; Kelly & Palumbi 2010). Indeed, coalescent analysis suggests that genetic structure observed in the more dispersive species results from ancient lineage divergence, which has been maintained with little gene flow between populations (McGovern et al. 2010). The role of gene flow, however, has not been explicitly tested from an oceanographic perspective. Here, we test the extent to which population genetic structure can be explained by recurrent and contemporary gene flow in this region, using a high-resolution oceanographic circulation model to simulate larval dispersal along the northwest coast of North America. We simulated gene flow of a marine invertebrate that has been the focus of numerous analyses of population genetic structure and coalescence using multiple genetic marker classes – the sea star Patiria miniata. This species has a long pelagic larval duration (6–10 weeks) with high potential for contemporary gene flow between populations, and provides a model

5038 J . M . S U N D A Y E T A L . for other long-duration passive dispersers. We first used gene flow simulations to test predictions of genetic structure against empirical estimates in P. miniata. We then used the model simulations to predict gene flow into and out of a previously unsampled location on the British Columbia coastline and tested these predictions with new genetic data for P. miniata from that location. We found that gene flow in the oceanographic model matched previously studied empirical patterns and provides relatively accurate predictions of population genetic structure in a previously unstudied part of the P. miniata geographic range.

Methods Study species We simulated population connectivity of the bat star, Patiria miniata, in which genetic population structure has been previously sampled using several classes of molecular marker: a mitochondrial sequence, seven nuclear microsatellites, two nuclear intron sequence markers (Keever et al. 2009), several anonymous nuclear sequence loci (McGovern et al. 2010) and sequences encoding gamete recognition genes (Sunday & Hart 2013). P. miniata adults live in the intertidal and shallow subtidal benthos, and their distribution is not continuous throughout their range. P. miniata are known to occur in California south of Fort Bragg, and on the outer coast of British Columbia, including the rocky shores of Haida Gwaii, and as far north as southeastern Alaska (Lambert 2000), but are absent from the coasts of Oregon and Washington. P. miniata adults spawn primarily in the summer (Rumrill 1989), followed by a planktotrophic larval phase estimated to be 6–10 weeks before settlement competency (Strathman 1987; Rumrill 1989; Basch 1996).

Oceanographic dispersal model and genetic model We used a buoyancy-driven, wind-driven, and tidally driven Lagrangian particle dispersal model (Blanton, 1995) for the coastal region of North America from California to southeastern Alaska to simulate larval dispersal to and from predefined coastal regions of equal size (Fig. 1). Ocean velocities were computed with a high-resolution diagnostic finite element model of the northeast Pacific Ocean forced with seasonal climatologies of temperature, salinity and wind stress (Foreman et al. 2008). We used the summer-averaged velocities to simulate the spawning season of P. miniata and defined a larval competency period as 30–70 days, the time within which we expected a larva to be capable of viable metamorphosis (Basch 1996). We held particles at

1 m depth from the surface, based on data for other echinoderm larvae (Pennington & Emlet 1986; Martin et al. 1997) and observations of negative-geotactic swimming behaviour in Pisaster ochraceus (Crawford & Jackson 2002) and Patiria miniata (JS personal observation). We focussed our simulations on the northern range of P. miniata, from southeastern Alaska to southern British Columbia, which has been most intensively genetically sampled. However, we also included a single region to represent the southern range of P.miniata in California. We used the most northerly region in which P. miniata have been observed in California, assuming it to be the likely point of recurrent connectivity with the northern range. From this simulation, we calculated the probability of gene flow to and from each region and used this connectivity matrix to inform a genetic model in which we simulated gene flow and genetic drift using multiple two-allele loci. After drift–migration equilibrium was reached (~500 generations), we calculated expected pairwise FST between each region as the mean across 50 independent runs, with each run representing an independent two-allele locus. We also calculated the dispersal distance of each simulated larva for comparison with empirically derived predictions of dispersal distance (Siegel et al. 2003). See the Appendix S1 (Supporting information) for a full description of the dispersal model, genetic model, sensitivity of pairwise FST to different model parameters and time, and calculation of particle dispersal distances.

Molecular sampling We calculated the mean pairwise FST between four regions of our model that were previously genetically sampled (Keever et al. 2009). We refer to these genetically sampled regions using bold capitals (BA, WH, HG, AK); five other regions that were included in the dispersal model but not in our population genetic sampling are identified using lower case italics (cvi, scc, shg, kit, pru; see Table 1 for definitions). We used microsatellite marker data for comparison with simulated data because of the closer match in allelic diversity between the seven microsatellite loci (3–12 alleles per locus) and allelic diversity in our connectivity model (two alleles per locus), in comparison with the much larger number of haplotypes observed within populations for mitochondrial and intron sequence data (Keever et al. 2009). If multiple sampled sites in Keever et al. (2009) were located within a single region in our model, we pooled these data prior to FST calculation. Specifically, the Bamfield and Tofino sites sampled in Keever et al. (2009) were pooled within the BA region, and the Sandspit, Moresby Camp and Taanuu sites were all pooled within the HG region. © 2014 John Wiley & Sons Ltd

O C E A N C I R C U L A T I O N M O D E L P R E D I C T S G E N E T I C S T R U C T U R E 5039 (a)

(b)

(c)

AK - Alaska pru - Prince Rupert kit - Kitimat HG - Haida Gwaii shg - Southern Haida Gwaii BB - Bella Bella scc - southern central coast WH - Winter Harbour cvi - Central Vancouver Island BA - Bamfield

(d) WA

OR

CA

Fig. 1 Dispersal of simulated larvae through an oceanographic circulation model on the west coast of British Columbia and southeastern Alaska. (a) Initial position of larvae (points) and larval recruitment areas (dark grey polygons) for each region. (b–d) Overlaid dispersal tracks of 1000 particles after 70 days. (c) Close-up within Queen Charlotte Sound and Hecate Strait. (d) Dispersal of particles from Fort Bragg, California.

Table 1 Connectivity matrix showing the proportion of particles moving from each dispersal area to each recruitment area, for pelagic larval competent to settle between 30–70 days, held at 1 m depth, using summer velocity fields To From Bamfield Central Vancouver Island Winter Harbour Southern Central Coast Bella Bella Southern Haida Gwaii Haida Gwaii Kitimat Prince Rupert Southeast Alaska

BA cvi WH scc BB shg HG kit pru AK

BA

cvi

WH

scc

BB

shg

HG

kit

pru

AK

0.755 0 0 0 0 0 0 0 0 0

0.155 0.4833 0.0917 0.1917 0 0 0 0 0 0

0 0.145 0.574 0.301 0 0 0 0 0 0

0 0 0.122 0.156 0.018 0 0 0 0 0

0 0 0 0 0.449 0 0 0.0259 0 0

0 0 0 0 0 0 0.855 0 0 0

0 0 0 0 0 0.534 0.0724 0 0 0

0 0 0 0 0.106 0.00097 0 0.427 0 0

0 0 0 0 0 0 0 0 0.6052 0.00099

0 0 0 0 0 0 0 0 0 0.037

Acronyms indicate locations described in left hand column and shown in Fig. 1.

New molecular samples We collected new microsatellite data for the same seven loci from a previously unsampled region, to test new predictions based on the dispersal model. Individuals were sampled from Bella Bella (BB, 52°100 N, 128°100 W), a location approximately midway between the northern and southern phylogeographic groups. We obtained tissue samples and microsatellite genotypes from 50 © 2014 John Wiley & Sons Ltd

individuals using the same methods as in Keever et al. (2009).

Analysis We tested for a correlative relationship between pairwise FST based on empirical genetic data (referred to as ‘observed FST’) and those predicted from the genetic model (referred to as ‘predicted FST’) using a Mantel test

5040 J . M . S U N D A Y E T A L . with 9999 permutations. We estimated the relationship using reduced major axis regression on raw data, using 1000 randomizations (Jensen et al. 2005). We also tested for a correlation between observed FST and the geographic distance between populations. Geographic distances were defined as the closest straight lines drawn between sample locations without crossing land. We ran Mantel tests and reduced major axis regressions first without the additional sample from Bella Bella to estimate the relationship predicted from existing data and then asked whether the pairwise FST values of the new samples fell along the predicted relationship. Reduced major axis regression was performed using the IBD Web Service (Jensen et al. 2005), and all other analyses were performed in R (R Development Core Team 2009). To specifically test the pattern of genetic structure between Bella Bella and previously sampled regions against the expectations from our model, it was necessary to propagate the uncertainty expected when using only seven loci and to scale model-generated predictions of pairwise FST to those calculated from microsatellite samples. We did this by iteratively running our genetics model through the connectivity matrix for seven independent model runs (simulating seven loci) and calculating pairwise FST between all populations and Bella Bella. We used the reduced major axis regression slope and intercept derived from previously sampled sites to convert model run FST (with two alleles per locus) to those expected among microsatellites (3–12 alleles per locus). We iterated this process 1000 times to generate a distribution of plausible FST estimates between Bella Bella and all other sampled populations. We used the clustering method in STRUCTURE (Pritchard et al. 2000) to heuristically search multilocus genotype groupings to determine how samples from Bella Bella were grouped with other populations. This method estimates the likely number of populations in the system (k) and assigns individuals a probability of belonging to each group, based on minimizing withingroup Hardy–Weinberg and linkage disequilibrium. Previous analyses (without Bella Bella) indicated two genotype clusters or groups located north and south of Queen Charlotte Sound (Keever et al. 2009). Our specific interest with the inclusion of Bella Bella samples was in differentiating among three hypotheses: clustering of Bella Bella individuals with northern samples; clustering of Bella Bella individuals with southern samples; or admixture of southern-like and northern-like genotypes in the Bella Bella sample. We used an admixture model without geographical priors, with allele frequencies correlated between populations, and admixture proportions estimated from the data and fixed for all populations. We ran Bayesian MCMC searches of

1 000 000 steps with a burn-in of 250 000. We carried out multiple (n = 7) model runs for each value of k up to k = 5 populations and used the method of Evanno et al. (2005) to find the best-fit value of k.

Results The dispersal model (using summer-averaged current velocities) predicted variable displacement of particles from their starting locations (Fig. 1). Within 70 days, particles from Vancouver Island tended to move northwards, although not as far as Haida Gwaii (Fig. 1b). Particles from the southern central coast of British Columbia (scc) moved distinctly offshore, towards Vancouver Island, and southwards along its west coast (Fig. 1c). Populations from Haida Gwaii were very isolated, with only a single particle (of 1000) moving from Kitimat to Haida Gwaii and no particles moving in the opposite direction (particle movement between kit and HG; Fig. 1c). These patterns are reflected in the connection probabilities of the connectivity matrix (Table 1). There was no successful particle movement across the geographic range disjunction between the northern (BC and Alaska) and the southern (California) parts of the P. miniata range. Although particles from California (representing the most northern known Californian population) moved almost due north (Fig. 1d), they did not travel far enough to reach any of the northern sites (Fig. 1d). Mean displacement distance across all populations based on a mean larval duration of 55 days was 64.2 km (Fig. 2). Displacement was notably higher for the California population, and the southern central coast population (Fig. 2). Only particles from the California population travelled distances as high as the scale predicted by the relationship in Siegel et al. (2003). The observed (microsatellite) and predicted (model) population pairwise FST values had very different magnitudes (Table 2), likely owing to the different number of alleles in the empirical microsatellite data (3–12) versus the simulated loci (two). Nevertheless, the matrices of observed and predicted population pairwise FST were significantly correlated in a Mantel test (r = 0.80, p = 0.042; Fig. 3a). Notably, these same pairwise FST values also showed a tight correlation with physical distance between sites, indicating isolation by distance, although this was marginally significant (r = 0.98, p = 0.081, Fig. 3b). The dispersal model also predicted FST reasonably well between the newly sampled Bella Bella region and other populations. Although the model overestimated connectivity between Bella Bella and Alaska, Haida Gwaii, and Winter Harbour, the predictions fell © 2014 John Wiley & Sons Ltd

O C E A N C I R C U L A T I O N M O D E L P R E D I C T S G E N E T I C S T R U C T U R E 5041 along the slope of the expected relationship (Fig. 3c); and three of four comparisons were within the error expected when sampling only seven loci (Fig. 4). The model predicted that the Bella Bella region would be genetically more similar to Winter Harbour (to the south) than to Haida Gwaii (to the north), despite a similar physical distance to each region (Figs 3c and 4). This prediction was based on the nonzero particle movement from Bella Bella to the adjacent southern central coast (BB to scc) and the relatively high particle movement from the southern central coast to Winter Harbour (BB to WH, Table 1). By contrast, the model predicted particle movement from Bella Bella north to Kitimat (BB to kit) but no dispersal from either Bella Bella or Kitimat to Haida Gwaii (HG or shg; Table 1). This prediction was borne out in the observed microsatellite frequency distributions, which revealed pairwise FST between Bella Bella and Haida Gwaii that was greater than that between Bella Bella and Winter Harbour (Fig. 3c).

Population clustering Based on the circulation model and connectivity matrix (Table 1), we predicted that propagules from Bella Bella would disperse to regions directly north and south (Kitimat and southern central coast) but from these regions would only spread south from the southern central coast to Winter Harbour. We therefore predicted that individuals in this region would be included in the Southern genotype cluster identified by Keever et al. (2009), including the Bamfield and Winter Harbour samples. Results from the STRUCTURE analysis indicated two clusters of microsatellite genotypes (Fig. 5), as in the previous STRUCTURE analyses of Keever et al. (2009). Most Bella Bella individuals were assigned to the group that occurs predominantly on Vancouver Island (Fig. 5). However, some Bella Bella individuals had high assignment probabilities to the second group that is characteristic of northern populations.

Discussion

Distance (km)

300 250 200 150 100 50 0

AK pru kit HG shg BB scc WH cvi BA Cali Population source Fig. 2 Density plots showing the distribution of absolute distances over which simulated larvae were displaced from each source site, midway through the competency period (larvae 55 days old). Dotted line shows mean distance travelled among simulated larvae, and dot-and-dash line shows predicted displacement for the same pelagic duration based on the relationship between pelagic larval duration and dispersal distances inferred from genetic data from 32 species within the range of these geographic locations (Siegel et al. 2003).

Our coupled oceanographic and genetic model produced predictions of population genetic structure that were well matched to observations. Although the distance between populations alone could explain the stepping-stone pattern of variation (strong but marginally significant correlation between pairwise FST and distance), our model critically demonstrates how this pattern can come to be, based solely on surface currents that retain long-lived larvae within nearby adjacent regions (100 k years, suggesting that these two groups were historically vicariant through several cycles of the Pleistocene glaciation (McGovern et al. 2010). This analysis estimated a low level of recurrent gene flow between the sampled populations in southeastern Alaska and Vancouver Island. Our model thus provides an explanation of how these two ancient lineages have come to be maintained despite the long pelagic duration for this species, via high self-recruitment and low, stepping-stone dispersal. Because we allowed our genetic model to run to drift–migration equilibrium, the source of the allelic variation that we simulated was arbitrary (we used two alleles of unknown genetic distance at 0.5 frequency in each population). We also held effective population sizes equal when this likely varies throughout the study region. Some of the discrepancy between our modelled and observed structure may therefore be attributable to the influence of ancient vicariance or to differences in effective population sizes in the natural populations. For example, the lasting influence of a north–south ancient vicariance may explain why our equilibrium model overpredicted FST between Bella Bella and regions to the north (Fig. 3). Overall, the model predicted relatively short distances of passive larval dispersal throughout the system (mean of 64.2 km) and high self-recruitment in most populations. We expect these to be overestimates of true dispersal distances because larvae can settle earlier than the mean duration used in the simulation (55 days), and because larval mortality would further bias realized dispersal in favour of earlier settlement after shorter dispersal distances. This finding is counter to general predictions for species with long pelagic durations (e.g. predictions from Siegel et al. 2003; see review in Cowen et al. 2000), but is consistent with the genetic pattern observed in P. miniata. The exception was the California-released propagules, which moved © 2014 John Wiley & Sons Ltd

distinctly northwards in our oceanographic model at higher relative velocities compared with those in British Columbia and Alaska. The short dispersal distances in British Columbia and southeastern Alaska suggest that the complex oceanographic features of this region may be distinctive in disrupting along-shore dispersal. Such oceanographic features may have a general role in the poor correspondence between pelagic larval duration and gene flow (Weersing & Toonen 2009; Selkoe & Toonen 2011) and the spatial pattern of genetic structure observed in other dispersive asteroids (Timmers et al. 2012). While the dispersal model had a reasonable ability to predict the genetic structure found in British Columbia and southeastern Alaska, it was not able to resolve the extremely high genetic similarity between southern British Columbia and California observed in Keever et al. (2009). Particles seeded in Vancouver Island did not travel far south, and although particles seeded in California travelled distinctly northwards with the greatest velocities, they did not move far enough to connect with the northern region of the P. miniata range. We suggest three alternative processes causing this discrepancy. First, connectivity may be greater in some extreme years, such as during El Nino Southern Oscillation (ENSO) events, which were not captured in our dispersal model. A northward coastal current is known to strengthen during ENSO years (Huyer et al. 2002) and has been invoked to explain sporadic and temporary occurrences of California-affiliated species on Vancouver Island (Schoener & Fluharty 1985; Behrens Yamada & Gillespie 2008). Indeed, a regional ocean circulation model predicted northwards dispersal from California 3.5 times greater during the late summer and fall in El Ni~ no years (See & Feist, 2010). The regularity of ENSO events lends support to this mechanism for the high similarity in P. miniata allele distributions among California and Vancouver Island for many unlinked loci (Keever et al. 2009). Second, this genetic pattern may be maintained by intermediate populations occurring

5044 J . M . S U N D A Y E T A L . along the Oregon or Washington coasts, which have either gone undetected, or have recently become extinct. Third, post-glacial colonization of Vancouver Island from California, together with large effective population sizes and slow genetic drift (Keever et al. 2009), may explain similarities of these two populations despite evidence of low contemporary gene flow between them. There are several simplifying assumptions inherent in the ocean circulation model. First, the circulation model used mean seasonal wind and buoyancy currents, without allowing fluctuations in velocities due to intra- or interannual variation (e.g. storms or ENSO events). Greater variance in dispersal distances could lead to greater overall connectivity. Second, the particle-tracking model did not include a stochastic element to allow diffusion of particles outside its deterministic track. This potentially limited the number of rare vagrants connecting certain pairs of populations, leading to an underestimate of connectivity. Third, we modelled dispersal of larvae at 1 m depth to simulate negatively geotactic larvae, but velocities will differ with depth, and dispersal of larvae of demersal or vertically migrating species would likely differ (Tankersley et al. 1995; Robinson et al. 2005). Nevertheless, the congruence between predictions generated from the model and observed genetic structure suggests that it provides reasonably robust predictions at the spatial scales we examined and warrants future use in predictive genetic studies in other taxa. One important application of our findings is how they inform a broader interpretation of the varying patterns of genetic structure recorded among marine fish and invertebrates in British Columbia and southeastern Alaska with prolonged planktonic development and high dispersal potential. For example, localized natural selection and active larval retention (directional swimming by large muscular larvae) have been suggested to explain the high genetic structure found in the rockfish, Sebastes helvomaculatus, with a genetic disjunction between Vancouver Island and Haida Gwaii similar to P. miniata (Rocha-Olivares & Vetter 1999). However, our results suggest that this structure may alternatively be caused by historical processes (such as vicariance) and maintained by oceanographic effects on passive dispersal alone (Hilbish 1996). By contrast, additional mechanisms are required to explain the lack of genetic structure in taxa with similar dispersal potentials in the same region, for example the sea star Pisaster ochraceus (Harley et al. 2006; Popovic et al. 2014), dispersing snails (Kyle & Boulding 2000; Marko 2004), chitons with dispersive larvae (Kelly & Eernisse 2007), intertidal copepods (Edmands 2001) and other dispersing species (Kelly & Palumbi 2010). The relatively low connectivity generated by our oceanographic model suggests that

the lack of genetic structure in other taxa is more likely attributable to historical processes such as relatively recent range expansions, or to large effective population sizes that limit the rate of lineage sorting through genetic drift. Coalescent methods that simultaneously estimate gene flow, the timing of isolation and effective population sizes could be used to independently test this hypothesis (Hart & Marko 2010). Our results have specific implications for marine spatial planning along the coastline of British Columbia. Canada is currently developing a bioregional network of marine protected areas (MPAs) to expand its protection of marine ecosystems. Ensuring demographic connectivity between MPAs was a key planning recommendation of a panel convened by the Canadian Parks and Wilderness Society (Jessen et al. 2011), echoing the broader consensus for the importance of MPA connectivity (Halpern 2003; UNEP-CBD 2008; Gaines et al. 2010; Hamilton et al. 2010). This panel recommended MPA distances of 20 to 200 km based on previous estimates of dispersal distances in other regions of the world (Jessen et al. 2011). However, our findings demonstrate that distances between MPAs should fall within the lower end of this range (20 km) to ensure connectivity between protected populations, even for species with the longest duration passive larval dispersal in the complex coastline of British Columbia (Fig. 2). In addition, some regions in our model had greater roles as ‘nodes’ within the dispersal network than others (Watson et al. 2011). Notably, the southern central coast region (scc) had high dispersal distances and was key in linking the central coast populations to those on the west coast of Vancouver Island. We can tentatively recommend the southern central coast region as a key connectivity node based on simulated dispersal patterns alone. However, these recommendations are based on simulations of summer-spawning surface dispersers, and exploring different seasons and depths of dispersal will be needed to generalize these recommendations to species with other life histories (e.g. Robinson et al. 2005). Finally, as sea surface temperature increases with climate warming and new habitats become available, modifications to this physically driven dispersal model may help to predict and understand patterns of new species colonization (Sorte 2013).

Acknowledgements We wish to thank Felix Breden, Michael Hellberg, Carson Keever, Susana Pati~ no, the Earth2Ocean lab, the Fab*lab and three anonymous reviewers for helpful feedback. Thanks to Morgan Hocking and Joel Harding for sample collection. JS, MH and WP were funded by the Natural Sciences and Engineering Research Council of Canada.

© 2014 John Wiley & Sons Ltd

O C E A N C I R C U L A T I O N M O D E L P R E D I C T S G E N E T I C S T R U C T U R E 5045

References Arndt A, Smith MJ (1998) Genetic diversity and population structure in two species of sea cucumber: differing patterns according to mode of development. Molecular Ecology, 7, 1053–1064. Basch LV (1996) Effects of algal and larval densities on development and survival of asteroid larvae. Marine Biology, 126, 693– 701. Baums IB, Paris CB, Cherubin LM (2006) A bio-oceanographic filter to larval dispersal in a reef-building coral. Limnology and Oceanography, 51, 1969–1981. Behrens Yamada S, Gillespie GE (2008) Will the European green crab (Carcinus maenas) persist in the Pacific Northwest? ICES Journal of Marine Science, 65, 725–729. Blanton BO (1995) DROG3D: User’s manual for 3-dimensional drogue tracking on a finite element grid with linear finite elements. Unpublished manuscript, Ocean Processes Numerical Methods Laboratory, University of North Carolina at Chapel Hill. Briggs JC (1974) Marine Zoogeography. McGraw-Hill, New York. Cowen RK, Lwiza KMM, Sponaugle S, Paris CB, Olson DB (2000) Connectivity of marine populations: open or closed? Science, 287, 857–859. Crawford BJ, Jackson D (2002) Effect of microgravity on the swimming behaviour of larvae of the starfish Pisaster ochraceus. Canadian Journal of Zoology/Revue Canadienne De Zoologie, 80, 2218–2225. Dodimead AJ, Favorite F, Hirano T (1963) Salmon of the North Pacific Ocean: Part II. Review of oceanography of the subarctic Pacific region: Int. N. Pacific Fisheries Commission, Vancouver B.C. 195 pp. Dupont L, Ellien C, Viard F (2007) Limits to gene flow in the slipper limpet Crepidula fornicata as revealed by microsatellite data and a larval dispersal model. Marine Ecology Progress Series, 349, 125–138. Edmands S (2001) Phylogeography of the intertidal copepod Tigriopus californicus reveals substantially reduced population differentiation at northern latitudes. Molecular Ecology, 10, 1743–1750. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 14, 2611–2620. Foreman MGG, Crawford WR, Cherniawsky JY, Galbraith J (2008) Dynamic ocean topography for the northeast Pacific and its continental margins. Geophysical Research Letters, 35, L22606, doi: 10.1029/2008GL035152. Foster NL, Paris CB, Kool JT et al. (2012) Connectivity of Caribbean coral populations: complementary insights from empirical and modelled gene flow. Molecular Ecology, 21, 1143–1157. Gaines SD, White C, Carr MH, Palumbi SR (2010) Designing marine reserve networks for both conservation and fisheries management. Proceedings of the National Academy of Sciences of the United States of America, 107, 18286–18293. Galindo HM, Olson DB, Palumbi SR (2006) Seascape genetics: a coupled oceanographic-genetic model predicts population structure of Caribbean corals. Current Biology, 16, 1622–1626. Galindo HM, Pfeiffer-Herbert AS, McManus MA, Chao Y, Chai F, Palumbi SR (2010) Seascape genetics along a steep cline: using genetic patterns to test predictions of marine larval dispersal. Molecular Ecology, 19, 3692–3707.

© 2014 John Wiley & Sons Ltd

Gilg MR, Hilbish TJ (2003) The geography of marine larval dispersal: coupling genetics with fine-scale physical oceanography. Ecology, 84, 2989–2998. Halpern BS (2003) The impact of marine reserves: do reserves work and does reserve size matter? Ecological Applications, 13, S117–S137. Hamilton SL, Caselle JE, Malone DP, Carr MH (2010) Incorporating biogeography into evaluations of the Channel Islands marine reserve network. Proceedings of the National Academy of Sciences of the United States of America, 107, 18272–18277. Harley CDG, Pankey MS, Wares JP, Grosberg RK, Wonham MJ (2006) Color polymorphism and genetic structure in the sea star Pisaster ochraceus. Biological Bulletin, 211, 248–262. Hart MW, Marko PB (2010) It’s about time: divergence, demography, and the evolution of developmental modes in marine invertebrates. Integrative and Comparative Biology, 50, 643–661. Hellberg ME (1996) Dependence of gene flow on geographic distance in two solitary corals with different larval dispersal capabilities. Evolution, 50, 1167–1175. Hey J, Nielsen R (2004) Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics, 167, 747–760. Hickerson MJ, Cunningham CW (2005) Contrasting quaternary histories in an ecologically divergent sister pair of low-dispersing intertidal fish (Xiphister) revealed by multilocus DNA analysis. Evolution, 59, 344–360. Hilbish TJ (1996) Population genetics of marine species: the interaction of natural selection and historically differentiated populations. Journal of Experimental Marine Biology and Ecology, 200, 67–83. Huyer A, Smith RL, Fleischbein J (2002) The coastal ocean off Oregon and northern California during the 1997–8 El Ni~ no. Progress in Oceanography, 54, 311–341. Irwin DE (2002) Phylogeographic breaks without geographic barriers to gene flow. Evolution, 56, 2383–2394. Jensen JL, Bohonak AJ, Kelley ST (2005) Isolation by distance, web service. BMC Genetics, 6, doi: 10.1186/1471-2156-6-13. Jessen S, Chan K, C^ ote I et al. (2011) Science-based Guidelines for MPAs and MPA Networks in Canada, pp. 58. Vancouver: Canadian Parks and Wilderness Society. Keever CC, Sunday J, Puritz JB et al. (2009) Discordant distribution of populations and genetic variation in a sea star with high dispersal potential. Evolution, 63, 3214–3227. Kelly RP, Eernisse DJ (2007) Southern hospitality: a latitudinal gradient in gene flow in the marine environment. Evolution, 61, 700–707. Kelly RP, Palumbi SR (2010) Genetic structure among 50 species of the northeastern Pacific rocky intertidal community. PLoS ONE, 5, e8594. Kenchington EL, Patwary MU, Zouros E, Bird CJ (2006) Genetic differentiation in relation to marine landscape in a broadcast-spawning bivalve mollusc (Placopecten magellanicus). Molecular Ecology, 15, 1781–1796. Kyle CJ, Boulding EG (2000) Comparative population genetic structure of marine gastropods (Littorina spp.) with and without pelagic larval dispersal. Marine Biology, 137, 835–845. Lambert P (2000) Sea Stars of British Columbia, Southeast Alaska, and Puget Sound. University of British Columbia Press, Vancouver, BC.

5046 J . M . S U N D A Y E T A L . Levin L (2006) Recent progress in understanding larval dispersal: new directions and digressions. Integrative and Comparative Biology, 46, 282–297. Marko PB (2004) ‘What’s larvae got to do with it?’ Disparate patterns of post-glacial population structure in two benthic marine gastropods with identical dispersal potential. Molecular Ecology, 13, 597–611. Marko PB, Hart MW (2011) The complex analytical landscape of gene flow inference. Trends in Ecology & Evolution, 26, 448–456. Martin D, Claret M, Pinedo S, Sarda R (1997) Vertical and spatial distribution of the near-shore littoral meroplankton off the Bay of Blanes (NW Mediterranean Sea). Journal of Plankton Research, 19, 2079–2089. McGovern TM, Keever CC, Saski CA, Hart MW, Marko PB (2010) Divergence genetics analysis reveals historical population genetic processes leading to contrasting phylogeographic patterns in co-distributed species. Molecular Ecology, 19, 5043–5060. Nielsen R, Wakeley J (2001) Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics, 158, 885–896. Palumbi SR (2004) Marine reserves and ocean neighborhoods: the spatial scale of marine populations and their management. Annual Review of Environment and Resources, 29, 31–68. Pennington JT, Emlet RB (1986) Ontogenetic and diel vertical migration of a planktonic echinoid larva, Dendraster excentricus (Eschscholtz) - occurrence, causes, and probably consequences. Journal of Experimental Marine Biology and Ecology, 104, 69–95. Pielou EC (1991) After the Ice Age: The Return of Life to Glaciated North America. Univ. Chicago Press, Chicago. Popovic I, Marko P, Wares JP, Hart MW (2014) Selection and demographic history shape the molecular evolution of the gamete compatibility protein bindin in Pisaster sea stars. Ecology and Evolution, 4, 1567–1588. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics, 155, 945–959. R Development Core Team (2009) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Robinson CLK, Morrison J, Foreman MGG (2005) Oceanographic connectivity among marine protected areas on the north coast of British Columbia, Canada. Canadian Journal of Fisheries and Aquatic Sciences, 62, 1350–1362. Rocha-Olivares A, Vetter RD (1999) Effects of oceanographic circulation on the gene flow, genetic structure, and phylogeography of the rosethorn rockfish (Sebastes helvomaculatus). Canadian Journal of Fisheries and Aquatic Sciences, 56, 803–813. Rumrill SS (1989) Population size structure, juvenile growth, and breeding periodicity of the sea star Asterina miniata in Barkley Sound, British Columbia. Marine Ecology-Progress Series, 56, 37–47. Schoener A, Fluharty DL (1985) Biological anomalies off Washington in 1982–83 and other major Ni~ no periods. In: El Ni~ no Effects in the Eastern Subarctic Pacific Ocean (eds Wooster WS, Fluharty EL), pp. 211–225. Wash. Sea Grant Program, Univ. of Washington, Seattle.

Schunter C, Carreras-Carbonell J, MacPherson E et al. (2011) Matching genetics with oceanography: directional gene flow in a Mediterranean fish species. Molecular Ecology, 20, 5167–5181. See KE, Feist BE (2010) Reconstructing the range expansion and subsequent invasion of introduced European green crab along the west coast of the United States. Biological Invasions, 12, 1305–1318. Selkoe KA, Toonen RJ (2011) Marine connectivity: a new look at pelagic larval duration and genetic metrics of dispersal. Marine Ecology Progress Series, 436, 291–305. Siegel D, Kinlan B, Gaylord B, Gaines S (2003) Lagrangian descriptions of marine larval dispersion. Marine Ecology Progress Series, 260, 83–96. Slatkin M (1987) Gene flow and the geographic structure of natural populations. Science, 236, 787–792. Sorte CJB (2013) Predicting persistence in a changing climate: flow direction and limitations to redistribution. Oikos, 122, 161–170. Strathman MF (1987) Reproduction and Development of Marine Invertebrates of the Northern Pacific Coast. Univ. Washington Press, Seattle. Sunday JM, Hart MW (2013) Sea star populations diverge by positive selection at a sperm-egg compatibility locus. Evolution and Ecology, 3, 640–654. Tankersley RA, McKelvey LM, Forward RB (1995) Responses of estuarine crab megalopae to pressure, salinity and light implications for flood tide transport. Marine Biology, 122, 391–400. Timmers MA, Bird CE, Skillings DJ, Smouse PE, Toonen RJ (2012) There’s no place like home: Crown-of-thorns outbreaks in the central pacific are regionally derived and independent events. PLoS ONE, 7, e31159. UNEP-CBD (2008) Decision adopted by the conference of the parties to the convention of biological diversity at its ninth meeting: (United Nations Environment Programme -Convention on Biological Diversity)UNEP/CBD/COP/DEC/IX/20. www.cbd.int/doc/decisions/cop-connect09/cop-09-dec-20-en.pdf. Watson JR, Siegel DA, Kendall BE, Mitarai S, Rassweiller A, Gaines SD (2011) Identifying critical regions in smallworld marine metapopulations. Proceedings of the National Academy of Sciences of the United States of America, 108, E907–E913. Weersing K, Toonen RJ (2009) Population genetics, larval dispersal, and connectivity in marine systems. Marine Ecology Progress Series, 393, 1–12. White C, Selkoe KA, Watson J, Siegel DA, Zacherl DC, Toonen RJ (2010) Ocean currents help explain population genetic structure. Proceedings of the Royal Society B-Biological Sciences, 277, 1685–1694. Whitlock MC, McCauley DE (1999) Indirect measures of gene flow and migration: F-ST not equal 1/(4Nm+1). Heredity, 82, 117–125. Wilson AB (2006) Genetic signature of recent glaciation on populations of a near-shore marine fish species (Syngnathus leptorhynchus). Molecular Ecology, 15, 1857–1871. Wright S (1969) Evolution and the Genetics of Populations. Variability Within and Between Populations. Univ. of Chicago Press, Chicago.

© 2014 John Wiley & Sons Ltd

O C E A N C I R C U L A T I O N M O D E L P R E D I C T S G E N E T I C S T R U C T U R E 5047

J.S. and M.H. designed research, I.P. provided new microsatellite data and analyses, W.P. provided geographic analyses, M.F. provided dispersal modelling tools, J.S. ran analyses and wrote the manuscript; all authors contributed to text and figures.

and introns from Keever et al. 2009, EF165733– EF165790, EF165792–EF165971, FJ939314–FJ939328, FJ850593–FJ850958, FJ850243–FJ850592; Bindin gene sequences from Sunday et al. 2013: KC503670– KC503757.

Supporting information Data accessibility All microsatellite data and locations of simulated larvae at 0, 40 and 70 days are archived on Dryad (doi:10.5061/dryad.t9j7 g). Genetic data from other relevant studies are available on GenBank: mtDNA

© 2014 John Wiley & Sons Ltd

Additional supporting information may be found in the online version of this article. Appendix S1 Supporting Methods. Fig. S1 Pairwise FST between two example population pairs with time, under different model parameters for population size and number of larvae successfully dispersing from each.