CORRECTED PROOF Predicting the dispersal of wild Pacific oysters Crassostrea gigas ( Thunberg , 1793 ) from an existing frontier population — a numerical study

Non-native populations of Pacific oysters Crassostrea gigas (Thunberg 1793) are established around the United Kingdom (UK), with two genetically different stocks originating from separate introductions to the UK and France. In this study, we use a coupled biophysical model to simulate Pacific oyster larval transport, in order to investigate the dispersal of the species from a known population near their northern limit on the west coast of the UK (in the Milford Haven Estuary). The model included a pelagic phase, simulating different swimming behaviours, and a settlement phase based on a hydrospatial substrate map. Following successful settlement elsewhere, subsequent releases simulated potential population spread over successive generations. Our results suggest that, should there be sufficiently warm sea temperatures to allow reproduction, dispersal away from Milford Haven Estuary would most be southeast ward towards the Bristol Channel; but dispersal north and west to Ireland is also possible, depending heavily on pelagic swimming behaviour. Seasonal modifications to circulation were less influential. Our study increases understanding of factors that contribute to oyster population spread, and suggests methods for improved management through numerical predictions.


Introduction
The introduction and establishment of non-native species, facilitated by the continued expansion of trade and transportation networks, has long been recognised as a global problem with important ecological and economic impacts (Grosholz 2002). In the marine environment, the transfer of aquaculture products is a significant vector for non-native species, second only to shipping in global importance (Molnar et al. 2008). Movement of the target aquaculture organism can lead directly to introduction and invasion of these species in the wild, although hitch-hiking species, whether on the cultured organism or on equipment, are an additional source of invaders.
The Pacific oyster Crassostrea gigas (Thunberg, 1793), native to East Asia, is an important aquaculture species that has been introduced to over 50 countries for commercial production. It was originally introduced into Europe from East Asia perhaps as early as the 16 th Century (O'Foighil et al. 1998;Mineur et al. 2014) and initial introductions led to naturalised populations in Portugal, with subsequent spread and establishment in France and other parts of Europe. Here, it is assumed that C. gigas is a synonym of the Portuguese oyster Crassostrea angulata (Lamark, 1819) (Batista et al. 2007), although phenotypic and genetic differences have been found (Batista et al. 2009;Zhong et al. 2014). In the early part of the 20 th century, C. gigas became the main cultivated species following the decline of the native flat oyster Ostrea edulis Linneus, 1758 (Mineur et al. 2014) and, throughout the 20 th century, C. gigas has been purposefully imported from different parts of the world to Europe to enhance aquaculture production (Miossec et al. 2009). As a consequence of importation of C. gigas, naturalised populations have established throughout Europe, most notably in France but also in many northern European countries, including Germany, Denmark, Sweden and Norway (Miossec et al. 2009, reviewed in Lallias et al. 2015. In the United Kingdom (UK), C. gigas was imported to Conwy, Wales, in 1965 for aquaculture. Although not expected to breed in the wild owing to low water temperatures, observations in 1989 and 1990 showed signs of natural spatfall (Spencer et al. 1994). There are now extensive natural beds of C. gigas in SE and SW England, with genetic work suggesting some could be as a result of natural dispersion of larvae from France (Child et al. 1995;Lallias et al. 2015).
The secondary spread of introduced non-natives, beyond the point of first introduction by human vectors, to the wider region is an integral part of the invasion process (Kolar and Lodge 2001;Sakai et al. 2001;Blackburn et al. 2011;Byers et al. 2015). A range of abiotic and biotic factors, such as water circulation or habitat availability, can act as barriers to dispersal and settlement of non-native species, which, together with species' life history traits, determine their range limits (Bohn et al. 2015). As a consequence of frequent introductions of C. gigas throughout Europe and subsequent establishment in the wild, there are now a number of isolated populations in which the limiting factors for secondary spread are unclear. Changes in environmental conditions as a consequence of changes in climate are predicted to enhance the likelihood of establishment and spread of non-native species following introduction (Stachowicz et al. 2002;Walther et al. 2009).
Currently, the northern limit of C. gigas in Europe is found at 59.2ºN in Espevik in Norway (Wrange et al. 2010). Elsewhere, northern boundaries are found at The Wash on the east coast of Great Britain (52.9ºN), Milford Haven Estuary on the southwest coast of Wales (51.7ºN), and Lough Swilly on the north coast of Ireland (55ºN) (Kochmann et al. 2013;Lallias et al. 2015). With the exception of the Lough Swilly population, it is not known whether these northern-most populations are self-sustaining or whether, because water temperature is too low, they are effectively sterile and seeded by larvae from distant populations. Large scale simulations by Thomas et al. (2016) using a combination of satellite derived sea surface temperature and a dynamic energy budget model address how environmental change over past decades has modified the extent of the reproductive niche of C. gigas in Europe. Using a sea surface temperature threshold for spawning of 18 °C, Thomas et al. (2016) demonstrated a 1,400-km northward shift in the boundary of C. gigas' ability to spawn over a 30-year period, up to 2003. This work shows the northern limit for spawning to occur as far north as the southern Norwegian coast, but only as far north as Cardigan Bay (just north of Milford Haven Estuary) on the west coast of Britain. Consideration of environmental control of secondary spread must also take into account global climate change.
Hydrodynamic models coupled to particle tracking models (PTMs) are routinely used to investigate larval dispersal. Using PTMs, large numbers of particles can be tracked in a Lagrangian framework to represent larval transport within a population. This framework allows the trajectories of individual particles to be recorded and later analysed. Brandt et al. (2008) used a coupled population abundance and particle tracking model to reproduce observations of C. gigas invasions in the German Wadden Sea. The simulations helped attribute the multi-year invasion pattern to a combination of external larval supply and larval hitch-hiking from local fishing activity, although with large model uncertainties and an assumption of passive larval dispersal. North et al. (2008) used models of Chesapeake Bay, USA, that did include larval swimming behaviour of two oyster species (Crassostrea virginica (Gmelin, 1791) and Crassostrea ariakensis (Fujita, 1913)) to show that different behaviours had a greater influence on dispersal and population connectivity than interannual differences in circulation.
Here, we consider the potential spread of C. gigas from an established population in Milford Haven Estuary (called Milford Haven hereafter), SW Wales, which is near to the northern limit of the expected reproductive niche (Thomas et al. 2016). By combining hydrodynamic simulations of the Irish Sea with a PTM which simulates the life cycle of C. gigas larvae, the potential dispersal of the population was predicted under different biophysical conditions. We are not aware of other studies in Europe that have used a similar approach. The virtual larvae were advected and mixed by the local hydrodynamics, and also transported due to their individual biological traits, controlled by sensory cues such as temperature, pressure, and light (Pedersen et al. 2003;Neill and Kaiser 2008;North et al. 2008). Because of uncertainty over larval behaviour in the water column, a sensitivity approach was taken in which different scenarios of larval vertical migration were modelled. The spread of a population occurs through a "stepping stone" process, whereby new populations are seeded from the source population, due to larval behaviour and larval transport along conduits of oceanographic currents . We therefore also examined the possibility of "second generation" spawning from these new populations, which would lead to a gradual spread, "third generation" and so on, around the coast (Kimura and Weiss 1964;Cowen et al. 2006).
We focused our study on offshore transport at shelf-scale, where we believe our model to be sufficiently robust to provide realistic simulations; it must be noted that coastal and estuarine larval transport simulations would require increased spatial and temporal model resolution, together with fine-scale data associated with boundary parameterisations. It is also worth noting that similar applications in deeper waters may require higher resolution in the vertical plane. The mesoscale Irish Sea contains tidallyenergetic regions with strong mixing (Xing and Davies 2001;Horsburgh and Hill 2003), as well as regions of seasonal stratification (Simpson and Hunter 1974;Hill et al. 1994) which generate density-driven currents that play an important role in larval transport over several weeks. Of particular importance for our simulations is the Celtic Sea front residual, which regulates particulate transport from Wales towards Ireland (Horsburgh et al. 1998;Robins et al. 2013). Because of the proximity of the Milford Haven to the front, it is expected to affect larval transport from populations within the estuary. In addition, the Bristol Channel, to the southeast of Milford Haven, has an extremely large tidal range (in excess of 12 m during spring tides) and asymmetrical tidal currents, which may affect the movement of larvae travelling within the channel.

Modelled species
Crassostrea gigas larvae are planktonic and move freely within the water column, having three stages with different durations (trochophore, veliger, pediveliger) before becoming competent to settle (Helm et al. 2004). The larval phase typically lasts 2-4 weeks, with the duration determined by physiological temperature. C. gigas is a euryhaline and eurythermic species, able to colonize any hard substrate, such as rocks, pebbles and bivalve shells (Kochmann et al. 2013). They settle in waters less than 40 m deep (FAO 2016), though more typically in intertidal and subtidal waters down to 15 m (Miossec et al. 2009;Nehring 2013), preferring sheltered areas such as estuaries and bays.

Hydrodynamic mode
The hydrodynamics of the Irish Sea were simulated using a three-dimensional (3D), free-surface, primitive equation, hydrostatic model (sbPOM; Jordi and Wang 2012). The Irish Sea model domain extended from 50.0ºN to 56.0ºN and from 9.0ºW to 2.7ºW ( Figure 1). Prognostic variables (e.g., velocity, surface elevation, temperature) were solved on an orthogonal, staggered Arakawa-C grid in the horizontal plane. The horizontal grid resolution is 1/60° longitude per 1/40° latitude, giving a mean cell size of approximately 1.85 × 2.2 km. The vertical plane was divided into 20 equally-segmented sigma layers giving a mean resolution of 4.3 m at mean sea level. The model incorporated a Level 2.5 turbulence closure scheme (Mellor and Yamada 1982).
A 7-month simulation was computed for the period 01 March to 30 September 1990. The model was forced at the open boundaries with the primary six tidal constituents, including the dominant semidiurnal M 2 (lunar) and S 2 (solar) constituents, but also N 2 , K 1 , O 1 , and P 1 ; generated from an outernested model (see Neill et al. 2010 for details of the outer model). The source of synoptic meteorological fields was the European Centre for Medium-Range Weather Forecasts-Interim reanalysis (Simmons et al. 2006), available at 3-h intervals at a spatial resolution of 1.5°. The initial month was computed in order for the density-driven currents to fully develop from initial temperature conditions (spatially constant at 8 °C), after which model output was generated every 30 minutes. Our simulation was considered a typical year in terms of decadal (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) variability of bed shear stress, used as a proxy for wind variability (Neill et al. 2010), and also Irish Sea temperatures ). Young and Holt (2007) reproduced Irish Sea temperatures and concluded that temperature variability is driven predominantly by local forcing (rather than boundary forcing), since the flushing time of the Irish Sea is typically 1-2 years, which is considerably longer than the response to seasonally varying meteorological forcing. Our Irish Sea modelled tidally-driven and density-driven circulation and temperatures have been quantitatively assessed against data from a variety of sources (see below and also Robins et al. 2013).
Comparing simulated elevations with 18 coastal stations within the domain, tidal analysis revealed that amplitudes and phases of the modelled M 2 and S 2 harmonics gave root-mean square errors of approximately 14 and 6 cm and 10° and 12°, respectively. A similar comparison for depth-averaged currents with 15 tidal current meters (see Davies and Jones 1990;Young et al. 2000) gave root-mean square errors of 0.05 m s -1 and 0.03 m s -1 in amplitude, and 15° and 10° in phase, respectively. Normalised root mean squared errors (scatter indices presented as percentages) for elevations and tidal velocities were less than 10%.
A comparison of the simulated sea surface temperatures at Milford Haven has been made with observed monthly mean temperatures at a nearby monitoring station (Angle Bay; CEFAS 2015), together with satellite data in that region (NASA 2015). The simulated temperatures were within recorded levels of natural variability over the period 1990 to 2007 (Table 1). Here, temperatures were representative of a colder-than-average period. Additionally, simulated temperatures were compared with monthly-averaged sea surface temperatures and temperatures at 37 m depth (measured at 54.5ºN, 4.5ºW (see Figure 1) and provided by Port Erin Marine Laboratory, University of Liverpool), giving root mean squared errors of less than 0.3 °C. Finally, simulated surface-bottom temperature differences and associated baroclinic residual flows were compared with another model (Xing and Davies 2001), showing good qualitative agreement of frontal positions, stratification strengths, and shelf-scale baroclinic flows (see Robins et al. 2013, Figure 2 and 3). Given the above baroclinic assessment of the model, we are confident that our spin-up period was adequate for the purpose of this study.

Particle Tracking Model
A particle tracking model (PTM) was developed to simulate the likely dispersal patterns of C. gigas larvae. The PTM simulates larvae as individual particles that are independently advected and mixed by the 3D velocity and turbulence fields produced by the hydrodynamic model. The PTM was run separately (off-line) so that the larvae could be assigned additional transport due to their biological behaviour over the course of the pelagic larval duration (PLD). The larvae were considered to be punctiform and their mass was assumed to have negligible influence relative to the water current speeds (i.e., no inertial effects of the larvae were included in the model). Larval advection was calculated by tri-linear interpolation of the simulated 3D velocity field and iterated forward in time using a PTM time step of 10 minutes.
The larvae are mixed locally through sub-grid-scale turbulent mixing, based on random displacement models (random walks), where the eastwards change in position, Δx (m), over a PTM time step ∆ (600 s) is given by (Proctor et al. 1994): where R is a random number in the range [0,1] and r is the standard deviation of Rcos(2πR) with a value of 1/√6. The expression (R/r)cos(2πR) thus has the necessary properties of a mean of zero and a standard deviation of unity (Ross and Sharples 2004). A similar expression for northwards diffusion is given by: Spatially and temporally varying horizontal diffusivities, K x and K y (m 2 s -1 ), were output from the hydrodynamic model and linearly interpolated in the same way as the velocities. A random displacement model (Visser 1997;North et al. 2006) was used to simulate sub-grid-scale vertical turbulence: where K z (m 2 s -1 ) is the vertical diffusivity, calculated by the hydrodynamic model and trilinearly interpolated to the particle position, z, in the PTM algorithm, and ⁄ is evaluated at depth z.

PTM simulations
Situated within Milford Haven, the population release site extended from 51.68ºN, 5.06ºW to 51.71ºN, 5.15ºW, which enclosed six model cells (Figure 1). For each simulation, 10 5 larvae were randomly distributed laterally and vertically within the release site. Larvae that were advected onto land were reflected back into the sea, so that no early settling/wastage occurred and larvae remained in circulation, maximising the statistical analysis of the results (see subsection Statistical analysis), as done by North et al. (2008) and Robins et al. (2013). For the same reason, no larval mortality was incorporated into the PTM, and the results were presented as a proportional analysis rather than predictions.
A series of spawning dates were selected in late July and August to coincide with the warmest seasonal temperatures: 21 July, 31 July, 10 August, 20 August, and 30 August. Particles were initially released at different times of the day to test sensitivity of spawning to the tidal cycle (i.e., during peak flood and peak ebb tides). However, dispersal patterns were similar in both cases. Similar results have been found in other studies, for example, by Kim et al. (2010) simulating Eastern Oyster larval transport, and by Coscia et al. (2013) modelling cockle larvae in the Milford Haven area. The results presented hereafter, (where particles were released at the start of each day), were representative of all spawning throughout the tidal cycle. Based on previous work (Coon et al. 1990;Collet et al. 1999;Troost 2010), a PLD of three weeks was chosen for larval development, with an additional stage (up to 7 days) during which the larvae were considered competent to settle. Each simulation consisted of three stages: Stage 1 (1.5 days), particles were fertilized gametes or trochophores that cannot swim and were transported passively by the water currents; Stage 2 (19.5 days), the larvae developed from veligers to pediveligers and may adopt vertical swimming behaviour; and Stage 3 (up to 7 days), the larvae were competent and ready to settle; they swim towards the bed and remain in suspension until suitable grounds to settle or become stranded after 7 days.
The movement of larvae in the vertical plane is an important factor in regulating larval dispersal . Unfortunately, there is insufficient information available for C. gigas to develop a specific swimming model during Stage 2. Consequently, a sensitivity approach was adopted, whereby three generalised swimming models were simulated for Stage 2: passive transport (where movement is controlled solely by the currents), tidal swimming (where larvae swim upwards during rising sea levels and downwards during falling sea levels), and diel swimming (where larvae swim downwards during the day and upwards during the night). These strategies also provided valuable insight into the variability of dispersal due to different behavioural patterns. As the larvae grow, they develop increasingly effective swimming capabilities (Troost et al. 2008). This was simulated by a linear increase in swimming speed during Stage 2, from zero (day 1.5) to 1-5×10 -3 m s -1 (day 21).
The simulated transport of larvae resulted in potential settlement only when the following conditions were met: (1) the simulation must have entered the settling period (Stage 3, days 21-28); (2) the larva must have reached the bottom model layer; (3) the substrate must be appropriate (e.g., rock and gravel; Coon et al. 1990); and (4) the water depth at that location must be less than 40 m (FAO 2016). Substrate type was based on an Irish Sea substrate map created using marine hydrospatial data (www.edina.ac.uk/digimap). If conditions 1-4 were satisfied, the probability of settlement was linearly increased from 5% (day 22) to 100% (day 28), since there is some evidence that C. gigas delay or improve their ability to settle as they age (Coon et al. 1990). However, C. gigas are commonly only found in the intertidal zone; therefore, our 40 m threshold based on the FAO (2016) fact sheet may overestimate successful settlement zones.
Finally, additional spawning sites were selected for second generation PTM simulations, where new oyster populations would most likely become established. The rationale for choosing the secondary generation sites was based on the simulated dispersal  patterns from the first generation simulations (from Milford Haven). Second generation releases were from regions where first generation settlement was > 5% and with close proximity to the coast (< 5 km). Only the behavioural traits that led to settlement were used for second generation spawning; for example, we simulated second generation spawning, with tidal migration only, from Watchlet, Devon, because the first generation dispersal of tidal larvae from Milford Haven successfully settled there (but not other behaviours). In total, six scenarios were simulated for second generation spawning (see Figure 1). From these six simulations, the procedure was repeated, resulting in one third generation simulation. The purpose of these simulations was to determine the possibility of population spread around the Irish Sea over successive years, via so-called "stepping-stone" spawning. Settlement areas were categorised into distinct geographical regions: Ireland, Fishguard, St. Brides Bay, Milford Haven area, Milford Haven, South Wales, North Bristol Channel, and South Bristol Channel (Figure 1). Connectivity was calculated as the proportion of all released larvae (10 5 ) which settled in each region at the end of each simulation.
In summary, we simulated three plausible behaviour scenarios for larval spawning from Milford Haven during 1990, each scenario consisting of 5 separate simulations of 10 5 larvae (Table 2). Additionally, stepping-stone simulations were performed for the 1990 conditions, totalling 22 PTM simulations.

Statistical analysis
Analysis of variance tests were performed to determine if dispersal distance and settlement success differ between sources of variation (i.e., larval behaviour and spawning season).
Univariate statistical analyses were calculated using the software package SPSS v14 and PRIMER 6.0 was used for multivariate statistics. An alpha value of 0.05 was used to determine statistical significance. Parametric ANOVA was used when the data satisfied the assumption of homogeneity of variance as determined through Levene's test; otherwise, the nonparametric Kruskal-Wallis test was used. Multivariate tests used were Multi-Dimensional Scaling (MDS) and Analysis of Similarities (ANOSIM) (Clarke and Gorley 2006).

Larval dispersal
We simulated variability in dispersal patterns of wild Pacific oyster larvae, spawned from Milford Haven (PTM- 1.1; Figures 2, 3). Variability was caused by simulating three different behaviour models (passive, tidal, diel), but also by changing the larval release dates (21 July to 30 August, 1990). Large tidal ranges in the region (4-6 m) produced strong tidal velocities around headlands (up to 2.5 m s -1 ), which quickly advected particles offshore from Milford Haven. Passive larvae ( Figure 2) were either advected north in the characteristic northwards tidal residuals, or westwards entrained in the Celtic Sea front residual current (see Horsburgh et al. 1998); both residuals being of the order 0.05-0.1 m s -1 . Smaller proportions of larvae were transported up to 50 km eastwards into the Bristol Channel, which largely displayed flood-tidedominant residuals. Tidal larvae often displayed bimodal dispersal, with large proportions of larvae being transported either north-westwards or eastwards ( Figure 2). As the Bristol Channel is flood-tidedominant, tidally-synchronised larvae can be transported large distances in the flood-tide direction.
Larval dispersal distances tended to be greater for passive and tidal behaviours (average radial dispersal away from release site of 25 km and 39 km, respectively), compared with the diel behaviour strategies (average radial dispersal of 13 km) (Figure 3). In fact, tidal populations consistently covered larger distances (population-averaged dispersal distances of 35 km) than the other behaviours. Univariate statistical analyses indicated that the differences in population-averaged dispersal distances between the three behaviour strategies were significant (ANOVA,  Table 2) and for each swimming strategy (passive, tidal, diel). Self-recruitment is denoted in the first column ("MH Estuary"). Table 3. Multivariate data analysis (ANOSIM), showing dissimilarity (Manhattan distance) between behaviour and settlement success. F 2,9 = 8.158, p = 0.006). There was no significant effect of season (i.e., release date) on dispersal distances (ANOVA, F 4,9 = 0.494, p = 0.74), although Figure 2 does suggest qualitative differences in dispersal among release dates.

Larval settlement and connectivity
Our simulations showed high self-recruitment at Milford Haven; passive, tidal, and diel strategies produced self-recruitment proportions of 64%, 61%, and 46%, respectively (Figure 4). Self-recruitment showed no significant differences between behaviour strategies (one-way ANOVA using release dates as replicates) although there was a noticeable difference in the proportion settling in the Milford Haven area ( Figure 4); 22% of diel swimming settled in this area compared with <2% in passive and tidal larvae. Larval wastage was between 18% and 28% depending on behaviour, either because larvae remained in deep water or only encountered unsuitable substrate. The ability of larvae derived from Milford Haven populations to settle far away from the South Wales area was considered unlikely, but possible. For example, simulations of passive larvae predicted relatively large proportions of settlement on the Irish coast (0.55% on average, but 2.73% on 31 July 1990), whereas other simulations predicted < 0.1% connectivity (tidal), or no connectivity (diel), with Ireland. Conversely, the tidal-behaviour model predicted 9% of larvae settled in the Bristol Channel (north and south), whereas other behaviour models predicted negligible settlement there. The geographic pattern of settlement among behaviour strategies was analysed using a multivariate approach, with release dates acting as replicates for each of the three behaviours and settlement percentages at each of the eight regions as multivariate response variables. ANOSIM showed a significant overall effect of larval behaviour strategy on settlement pattern (R = 0.18, p = 0.029); pair-wise tests revealed a statistically significant dissimilarity only between the tidal and diel behaviours (Table 3).

Population spread
New oyster populations could potentially become established in the following regions identified by successful settlement from PTM-1.1: St. Govan's Head, Caldey Island, Ramsey Island, and Watchet ( Figure 1). We performed second-generation simulations (PTM-2.1-2.3; Table 2) from these locations, but  Table 2). only for the behaviour strategies that led to their seeding. One of the most significant differences of the second-generation releases was increased larval wastage -60% on average (Figure 5a), compared with 25% for PTM-1.1. This is because the secondgeneration sites are more exposed to strong tidal currents, being located near headlands or open coasts. Consequently, population-averaged self-recruitment was reduced from 57% (PTM-1.1; Figure 4) to 15% ( Figure 5a). However, proportions of connectivity increased for all second-generation releases, except from Watchet, since the tidally-synchronised larvae from Watchet travelled east up the Bristol Channel, which is mainly unsuitable for larval settlement being comprised of soft sediment. In particular, significant variance in connectivity between the established Milford Haven population and the second-generation populations were calculated for passive strategies using the nonparametric Kruskal-Wallis test (South Wales populations: χ 2 = 13.491, p = 0.004; Ramsey Island population: χ 2 = 3.962, p = 0.047). The Irish coast was reached by Ramsey Island passive releases, in greater proportions (4% on average, but up to 10.2% connectivity) than from Milford Haven (up to 2.7% connectivity with Ireland). In all secondary simulations, a small proportion (average of 1%) of larvae settled back in Milford Haven (Figure 5a). Second-generation dispersal distances increased significantly for passive larvae, from 23.7 km to 46.7 km (ANOVA, F 3,9 = 7.533, p = 0.002), but not for the other strategies (results not shown). Again, this can be attributed to exposure to strong tidal currents.
From the six second-generation simulations (Figure 5a), larvae released from both St. Govan's Head and Caldey Island led to high settlement in our defined "South Wales" coastal region. Specifically, high settlement of passive larvae occurred in Loughor Estuary, which matched our criteria for steppingstone migration (settlement > 5% and distance to coast < 5 km). Note: the high settlement from St. Governs-diel and Caldey-passive larvae (Figure 6a) did not settle in the Loughor Estuary in high enough proportions for our criteria, but did settle elsewhere along the coast. Therefore, we chose a third-generation simulation from the Loughor Estuary (PTM-3.1; Table 2). This passive simulation showed high retention (36%) and connectivity (54%), compared with previous simulations, and ~10% larval wastage (Figure 5b). The predominant larval trajectory was southeast, with high connectivity in the Bristol Channel and near Swansea. This scenario sustains the trend of dispersal along the coast of South Wales under the passive strategy. The Loughor being a sheltered estuary, dispersal distances were overall shorter than from (exposed) second-generation sites.

Discussion
Our results suggest that current northwards range expansion of C. gigas on the western UK coast is constrained by both physiological (i.e., swimming behaviour) and dispersal barriers. Successful reproduction of Milford Haven populations is more likely to lead to population expansion to the southeast, and also across to southern Ireland. Although few larvae actually reached the Irish coast, being stranded instead offshore, our results suggest that under favourable conditions (e.g., easterly winds), larvae are perhaps capable of settling on the Irish coast in significant numbers. Northwards range expansion would be plausible if northerly residual currents were to persist in the future, which is likely to be the case, based on simulations by Olbert et al. (2012). However, one of the main climate change stressors that may affect larval densities, growth, and mortality is reported to be ocean acidification; although impacts are thought to be species-specific (e.g., Miller et al. 2009) and species adaptation is presently unknown. Because of lack of projected data on ocean acidification, we did not take this into account in our study.

Temperature thresholds for larval development
The development time for larvae can be defined as the cumulative degree days above a minimum threshold temperature for gametogenesis: , where is the number of days required for to attain a ripe state, is the ambient temperature and is a threshold temperature. For spawning of Pacific oysters to occur, approximately 600 cumulative degree days above a threshold temperature of 10.55 °C are required (Mann 1979;Syvret et al. 2008). Based on these assumptions, Table 1 suggests that Pacific oysters will not successfully develop and spawn at the northern boundary of the species in western UK. Further, if spawning was to occur as per our simulations, our model calculates that insufficient cumulative degree days would be accumulated by the larvae for metamorphosis and settlement (based on Mann's (1979) value for of 10.55 °C), irrespective of the direction of dispersal. It is therefore unlikely that wild oysters in Milford Haven are the offspring of farmed oysters at that site, if the above assumptions are correct. This conclusion is supported by Lallias et al. (2015), who found that abandoned farmed oysters at Milford Haven were from a different genetic cluster to a wild population collected just a few metres apart. Although there is no firm evidence, Lallias et al. (2015) suggest that the origin of the wild population could be from unregulated or accidental introductions, or from larval dispersal from Brittany. It is also worth noting that the size range of samples taken from the Milford Haven area (N. Zwerschke, Université de Nantes, Nantes, unpublished data) indicate that spatfall had not occurred for some years, being either a sporadic or a one-off event.
Pacific oysters have become established in areas where temperatures regimes were thought to be unfavourable to their life cycle (Wehrmann et al. 2000). Therefore, before concluding that the wild population found in Milford Haven is functionally sterile, it is worth revisiting the assumptions about temperatures. Fabioux et al. (2005) found that oysters could reach maturity when conditioned at temperatures as low as 8 °C, suggesting that Mann's (1979) threshold temperature for gametogenesis of 10.55 °C is too high. Chávez-Villalba et al. (2002) found that, in France, northerly populations matured more quickly than those from the south when they were conditioned in a common environment, suggesting genetic differences may exist. The number of degree days required for maturation appears to vary between studies, although it is difficult to determine whether conditioning may have commenced prior to some trials (Syvret et al. 2008). A number of papers report that spawning at a threshold temperature of 16 °C, whilst often considered suboptimal, is possible (e.g., Gillespie et al. 2012;Norgard et al. 2014). Following spawning, temperature requirements for successful metamorphosis vary, with Syvret et al. (2008) suggesting that a further 225 degree days above a threshold of 10.55 °C are required, and other studies reporting that 18 °C for two weeks or 20 °C for three weeks are required-equivalent to 104 and 199 degree days, respectively, using Mann's threshold value for gametogenesis (Gillespie et al. 2012). Whilst larvae are reported not to survive well below 15 °C (Gillespie et al. 2012), they have been observed to survive to metamorphosis and settlement at 14 °C in laboratory conditions (P. Boudry, IFREMER, Brest, pers. comm.). In the wild, larvae have been recorded at low densities at temperatures as low as 13 °C around Vladivostok, in their native range, although it is not clear whether these early spawned larvae go on to successfully settle as spat (Kulikova et al. 2015). However, in the wild, the increased larval period would make them more vulnerable to predators (Pauley et al. 1988). Inter-individual variation in degree days required for maturation is observed in hatcheries as not all larvae metamorphose or settle at the same time. Quality of nutrition can also influence the duration of the larval phase (Rico-Villa et al. 2006). Once competent, Coon et al. (1990) found that in the laboratory larvae were able to delay metamorphosis and remain competent for at least a month at 23 °C. However, it is likely that larvae in the wild would be subject to settlement or metamorphosis cues within that time.
The likelihood of successful spawning at the northerly limit of C. gigas in the Irish Sea may increase in the future, as sea temperatures rise. Diederich et al. (2005), in their study of the establishment of C. gigas in the Wadden Sea, observed that successful spatfall and settlement occurred in years when temperatures were higher than average during July and August. Olbert et al. (2012) predicted that Irish Sea surface temperatures will increase by 1.9 °C during the 21 st century, based on projected global warming and sea-level rise of 0.47 m. If sea temperatures were to increase, based on continuing these projections, then regular spatfall could be achieved. As well as there being sufficient degree days in this scenario, temperatures would be high enough for larvae to reach spatfall with good survival rates.

Dispersal and settlement
Larval dispersal is controlled by hydrodynamics and assumed larval behaviour. There is varied evidence for the existence of vertical migration behaviour in oyster species, such as C. virginica (e.g., North et al. 2008) but sufficient doubt to justify examining different migration scenarios. In our study, passive larvae dispersed to Southern Irish Sea coasts. The tidally-synchronised swimming pattern predicted that new populations could establish themselves in the Bristol Channel. Larvae with diel migration remained within the Milford Haven area. Of course, variability existed within this rather simplistic dispersal description; for example, due to natural seasonal variability in the hydrodynamics, such as the development of the Celtic Sea front during warmer months. In tidally energetic regions like the Irish Sea, the interaction between strategy and oceanography is most prominent; for example, regions of tidal asymmetry (often coastal inlets) will greatly influence tidally-synchronised larval dispersal .
Based on our present-day biophysical modelling, larvae potentially released from Milford Haven were most likely to remain in the South Wales area, either within the estuary itself, or in the surrounding waters. No larvae settled in the English Channel; therefore, the Milford Haven oyster population is unlikely to be a source population for its nearest neighbouring populations to the south (i.e., Cornwall). Whilst a small number of larvae settled to the north of Milford Haven, most northwards dispersal ended in wastage. This is consistent with observations on other non-native species such as Austrominius modestus (Darwin, 1854) which, failed to spread northwards from Milford Haven when it first arrived (Crisp 1958); and Crepidula fornicata (Linneus, 1758), which forms an extensive breeding population in Milford Haven, but has not extended substantially northward over a 50 year period (Bohn et al. 2012).
Far-field connectivity, for example between South Wales and Ireland, has been shown to be possible under some scenarios using a combination of modelling and genetic evidence in the cockle Cerastoderma edule (Linneus, 1758) (Coscia et al. 2013). We simulated low, but potentially biologically important, proportions of oyster larvae settling on the Irish coast for passive (2.73%) and tidal (<0.1%) larva. The Pacific oyster is a highly fecund species, with a single female capable of producing 50-100 million eggs per spawning (Quayle 1969). Thus infrequent colonisation events from across the Irish Sea could lead to widespread establishment. The lack of far-field connectivity of diel larvae (no diel larvae reached Ireland) may be attributed to desynchronisation between the swimming behaviour and the tidal stream currents, in addition to large amounts of time spent in weak bottom waters.

Population spread
The ability of populations of sessile invertebrate, such as Pacific oysters, to spread over large areas is conditioned by the establishment of new populations at a distance, where reproductive settlement success replenishes mortality (Kimura and Weiss 1964). We found a general pattern in which initial spread from Milford Haven was limited owing to high levels of self-recruitment, but simulation of larval release from newly colonised sites led to greater connectivity with sites further afield, most likely owing to greater tidal energy in these secondary spawning sites. The model therefore revealed a trade-off between exposed sites that have higher connectivity and sheltered sites having higher self-recruitment. Connectivity from secondary sites back to Milford Haven is unlikely; the original population being almost entirely seeded by self-recruitment. Although for our model most northerly dispersal was simulated as larval wastage, the results show that it is physically plausible that northwards migration could eventually occur into Cardigan Bay (which has suitable substrates), provided temperature regimes would allow oysters to spawn and disperse from there.
Our multi-generation spread "trail" simulations ended with a third generation site (Loughor Estuary). The high level of self-recruitment, its sheltered location, and potentially higher temperatures (compared with Milford Haven) would make it a favourable location for the development of a new oyster population. More importantly, the directions of dispersal and magnitude of settlement support the view that oysters will gradually spread eastwards around the coasts of South Wales. Over three generations, certain dispersal paths, with Milford Haven as an initial population in 1990 and under the different swimming strategies, seem likely ( Figure 6). The existing C. gigas population at Milford Haven, if sea temperature was warm enough for reproduction to occur, could disperse around the south coast of Wales and the Bristol Channel, with the possibility of new wild populations establishing in South Wales, Cardigan Bay, and on the north coasts of Devon and Somerset.

Future recommendations
Future simulations should account for the inter-annual variability in the atmospheric climate. For instance, the ECMWF-ERA-Interim reanalysis showed that, for the decade 1989-1998, there was considerable variation in the winter-mean (December-March) NAO index, ranging from −2.32 (anticyclonic, cold and dry conditions) to +2.44 (strong westerlies with warm and wet conditions) (Neill et al. 2010). However, Neill et al. (2010) stated that although variability was high over this decade, annual residual bed shear stress generally remained similar, regardless of extreme wind conditions. Young and Holt (2007) showed that, over the period 1960-1999, there has been a long-term warming trend in the Irish Sea (+0.018 °C yr -1 ). Both mean-summer and meanwinter temperatures varied by 1 °C, with the onset and breakdown of stratification variable by up to one month, potentially leading to variability in densitydriven flows and larval trajectories. Simulations that include the above variabilities, together with downscaled climate change projections, and performed on a high resolution grid (e.g., 1 km 2 cells), would greatly improve our knowledge of the variability of larval dispersal. However, a lack of knowledge of larval behaviour for many species is currently the most limiting factor to improved dispersal simulations.

Conclusion
We present a modelling approach (particle tracking) to simulate larval transport of a non-native population of Pacific oyster, at their northern limit on the west coast of the UK. Presently, it is not clear whether wild populations have spread from other wild populations, or whether they are seeded from nearby oyster farms. We have investigated how wild populations can spread via the stepping-stone dispersal mechanism. We show significant variability in predicted dispersal patterns, based on a range of different larval swimming behaviours-highlighting the importance of improving knowledge of larval behaviour for accurate predictions of population spread. Seasonal modifications to circulation were less influential, although we did not investigate interannual variations or climate change. Nevertheless, we present a method for predicting the likely pathways of population spread over successive generations. Assuming adequate temperature condition for spawning and larval development, our work shows population spread from Milford Haven would most likely occur southeastward along the Welsh coast towards the Bristol Channel, but dispersal north and west to Ireland is also possible. Since environmental agencies are concerned about population spread in all directions, and regulators want to know whether sterile oysters should be farmed in new sites, this work demonstrates that modelling approaches such as ours could help with the decision process.