Spill-over from aquaculture may provide a larval subsidy for the restoration of mussel reefs

Worldwide bivalve aquaculture is expanding rapidly. Simultaneously, there has been a loss of natural bivalve reefs due to anthropogenic activities. As bivalve reefs support several ecosystem functions disproportionate to the area they cover, there is interest in their restoration. The Firth of Thames (FoT) in northern New Zealand once supported dense populations of green lipped musselsPerna canaliculus, which were extirpated by a dredge fishery in the mid-20thcentury. Efforts to restore these biogenic habitats are underway. The largest standing populations of this species in the area currently exist in aquaculture. This study aimed to determine if larval spill-over from aquaculture can provide a larval subsidy to bivalve reef restoration efforts in the FoT. We used a combination of trace elemental fingerprinting and biophysical modelling techniques to determine patterns of larval dispersal in the area. Results of both approaches indicated that the larval pool in the area is well mixed with larvae produced at aquaculture locations capable of settling throughout the study area. Overall this shows, for the first time, that larval spill-over from aquaculture may provide a subsidy to restoration efforts and assist with establishing sustainable populations. When determining restoration locations, the potential for aquaculture populations to act as a larval source should be explicitly considered. Conversely, when considering the location of new aquaculture sites, the consequences of larval spill-over to surrounding wild populations should be assessed. We recommend that restoration efforts and aquaculture be carefully integrated in a network approach which could provide both ecological and economic benefits.


INTRODUCTION
At a global scale, coastal shellfish aquaculture is rapidly expanding in order to help meet the protein needs of a growing human population (FAO 2018). As cultured biomass of shellfish increases, it will inevitably affect ecosystem dynamics in the local area in which this culturing is taking place (van der Schatte ABSTRACT: Worldwide bivalve aquaculture is expanding rapidly. Simultaneously, there has been a loss of natural bivalve reefs due to anthropogenic activities. As bivalve reefs support several ecosystem functions disproportionate to the area they cover, there is interest in their restoration. The Firth of Thames (FoT) in northern New Zealand once supported dense populations of green lipped mussels Perna canaliculus, which were extirpated by a dredge fishery in the mid-20 th century. Efforts to restore these biogenic habitats are underway. The largest standing populations of this species in the area currently exist in aquaculture. This study aimed to determine if larval spill-over from aquaculture can provide a larval subsidy to bivalve reef restoration efforts in the FoT. We used a combination of trace elemental fingerprinting and biophysical modelling techniques to determine patterns of larval dispersal in the area. Results of both approaches indicated that the larval pool in the area is well mixed with larvae produced at aquaculture locations capable of settling throughout the study area. Overall this shows, for the first time, that larval spill-over from aquaculture may provide a subsidy to restoration efforts and assist with establishing sustainable populations. When determining restoration locations, the potential for aquaculture populations to act as a larval source should be explicitly considered. Conversely, when considering the location of new aquaculture sites, the consequences of larval spill-over to surrounding wild populations should be assessed. We recommend that restoration efforts and aquaculture be carefully integrated in a network approach which could provide both ecological and economic benefits. alteration of flow regimes by the structures on which aquaculture is taking place (Gibbs et al. 1991, Plew 2011. Effects perceived as positive may involve nutrient remediation in eutrophied systems (Petersen et al. 2014, Nielsen et al. 2016 or providing habitat for other species (McLeod et al. 2014, Fariñas-Franco & Roberts 2018. Although often overlooked in studies investigating the effects of aquaculture on the local environment (see reviews by Newell 2004, van der Schatte Olivier et al. 2020, larval production from these dense aggregations of bivalves may contribute to regional larval pools and thereby modify population dynamics in the area in which they are growing. Knowledge of potential larval dispersal pathways from these cultured populations is essential to gain an understanding of how larvae from aquaculture may interact with local wild bivalve populations. Understanding these interactions is of particular interest in areas where natural populations of bivalves have been depleted but cultured populations exist; e.g. the Firth of Thames (FoT) in northern New Zealand for green lipped mussel or Chesapeake Bay for Virginia oysters (Beck et al. 2011, Turner et al. 2019.
Shellfish beds perform a number of ecosystem functions disproportionate to their size (Lundquist et al. 2017). They contribute to water filtration (zu Ermgassen et al. 2013, Ehrich & Harris 2015, nutrient cycling (Giles & Pilditch 2006, Nielsen et al. 2016, provide habitat structure (Dumbauld et al. 2009), and can influence food web dynamics and biodiversity (McLeod et al. 2014, Fariñas-Franco & Roberts 2018. However, several common anthropogenic stressors are recognised worldwide which have resulted in degradation or removal of these biogenic habitats with a corresponding loss of ecosystem functions (Beck et al. 2011, Alleway & Connell 2015. These stressors include sedimentation, eutrophication, overfishing, and dredging. Accordingly, there is interest in restoring degraded populations and the ecosystem functions they provide (e.g. French McCay et al. 2003, Coen et al. 2007, North et al. 2010. In order for restoration programmes to have the highest chances of success, an understanding of larval dispersal and population connectivity is essential (Lipcius et al. 2008). Understanding larval dispersal dynamics allows us to determine the selfsustainability of restored populations and also select optimal restoration locations (Elsäßer et al. 2013, Tettelbach et al. 2013.
Although it is important to track the dispersal of larvae, in situ tracking of larval movement is complex due to the small size, high mortality, and large scales over which larval transport is possible (Gawarkiewicz et al. 2007, Cook et al. 2014. Thus, several techniques have been developed which allow larval dispersal and population connectivity to be indirectly inferred (reviewed by Kool et al. 2013). Currently, biophysical models (e.g. Cowen et al. 2000, Werner et al. 2007, Elsäßer et al. 2013, Thomas et al. 2016 and trace elemental fingerprinting (e.g. Thorrold et al. 2002, Becker et al. 2007, Ricardo et al. 2015, Gomes et al. 2016) are among the most widely applied techniques used to infer the dispersal of larvae. Biophysical modelling couples oceanographic data (e.g. ocean currents and temperature) with information on the biology and behaviour of larvae to estimate larval dispersal over spatial and temporal scales that would be impossible to measure empirically. Biophysical models have been used in a number of management scenarios, such as the design of marine reserves (Puckett et al. 2014), examining larval export from marine reserves (Le Port et al. 2014), determining restoration locations (Tettelbach et al. 2013, predicting the spread of invasive species (Inglis et al. 2006), and projecting changes under future climate scenarios (Cetina-Heredia et al. 2015, Coleman et al. 2017. Trace elemental fingerprinting is based on the premise that carbonate structures such as shell or otolith deposited in water masses with differing physiochemical properties reflect these differences in the trace elemental composition of these structures (e.g. Cathey et al. 2012, Kroll et al. 2016. Through the sequential analysis of these structures, the movement of individuals can be inferred (e.g. Lazareth et al. 2003, Kroll et al. 2018. Each method of indirectly estimating larval dispersal, however, comes with intrinsic uncertainties (Ashford et al. 2010, Nolasco et al. 2018. Biophysical particle tracking models, for example, are limited by the availability and accuracy of the hydrodynamic and biological data used in their parameterisation (Cetina-Heredia et al. 2019). In contrast, trace elemental fingerprinting may be limited by the characterisation of all potential natal locations (Werner et al. 2007). Therefore, the incorporation of larval dispersal estimates obtained though different methods is likely to provide a more accurate picture of larval dispersal. Here, we incorporated results from both trace elemental fingerprinting and biophysical models to track the dispersal of green lipped mussels Perna canaliculus (Gmelin, 1791), hereafter mussels, in the FoT in northern New Zealand, an area with an active shellfish restoration effort (Wilcox et al. 2018).
The FoT once supported dense populations of mussels, which were estimated to cover up to 1300 km 2 ; these were almost completely extirpated by a dredge fishery in the 1960s (Paul 2012). Despite the cessation of fishing 50 yr ago, these populations have not recovered and only a few small remnant wild mussel beds remain, which cover an area of approximately 0.64 km 2 (Morrison 2002, McLeod 2009). Currently, the largest populations of mussels in the FoT are those farmed on longline aquaculture, which covers an area of approximately 29 km 2 and produced 27 196 t of harvested mussels in 2016 (Hauraki Gulf Forum 2017). There are currently active restoration efforts in the area which focus on the translocation of adults and juveniles, primarily from aquaculture, with the aim of establishing self-sufficient populations and recovering lost ecosystem function (Wilcox et al. 2018). Understanding the potential for larval spill-over from aquaculture to contribute to restoration efforts will provide important information for the design of restoration programmes in areas where aquaculture is present, as it will in form locations for mussel translocations.
The aim of this study was to examine the potential for larval spill-over of mussels from aquaculture to provide a larval subsidy to the restoration efforts in the FoT. This was achieved using a combination of elemental fingerprinting techniques and biophysical modelling approaches. We addressed the following specific questions: (1) What is the potential for dispersal of mussles from aquaculture locations throughout the FoT? (2) How do larval dispersal patterns generated by 2 different methods, biophysical modelling and trace elemental fingerprinting, differ? (3) What are the overall implications of any larval spill-over for restoration and aquaculture?

Study species
Green lipped mussels Perna canaliculus are found in a variety of habitats throughout New Zealand, both subtidally and intertidally (Morton & Miller 1973). Mussels are broadcast spawners with peak spawning occurring from late in the austral spring to early autumn (Alfaro et al. 2001), have a pelagic larval duration of between 3 and 5 wk (Hayden 1995), and prefer to settle on filamentous substrates such as macroalgae and hydroids (Alfaro & Jeffs 2002, Alfaro et al. 2004. Settlement is complete once byssal threads have attached, after which mussels are known as plantigrades or spat. Secondary settlement may then occur through bysso-pelagic drifting from primary settlement substrate into adult beds (Jeffs et al. 1999). The distances over which secondary dispersal takes place, however, is likely to be orders of magnitude lower than primary dispersal in the plankton (Lane et al. 1985, Le Corre et al. 2013, Pilditch et al. 2015.

Study area
The FoT is a large estuarine embayment in north eastern New Zealand, approximately 30 km long and 20 km wide (Fig. 1). Circulation in the region is tidally driven and predominantly flows north−south, with current speeds on the ebb and flood tides generally well balanced (Supplement 1 at www.int-res. com/ articles/ suppl/ q012 p231 _ supp.pdf [all supplements], Black et al. 2000). East−west flows also occur but are typically around 10 times weaker than the north−south currents (Oldman et al. 2007). The FoT is generally thought to be well mixed with limited seasonal stratification in late summer (Black et al. 2000, Zeldis et al. 2005. Primary freshwater inputs are from the Waihou and Piako Rivers, which drain into the south eastern FoT. These rivers have a combined mean flow of 89 m 3 s −1 (O'Callaghan & Stevens 2017). The area zoned for mussel aquaculture occupies an area of 29 km 2 and accounts for 30% of New Zealand's Greenshell TM (the trade name of P. canaliculus) mussel production (Hauraki Gulf Forum 2017). These farms are primarily located on the northern end of the Coromandel Peninsula centred on Wilson Bay. There are additional farms located within and around Coromandel Harbour and Manaia Harbour and on the western FoT at Waimang Point.

Biophysical particle tracking
2.3.1. Hydrodynamic model An 11 yr (1995An 11 yr ( −2005 hindcast hydrodynamic model configuration of the regional ocean modelling system (ROMS), developed by the Meteorological Service of New Zealand (MetOcean Division), was used in this study. For a full description of this model and its validation, see Supplement 1. Briefly, the model has a spatial grid resolution of 250 × 250 m over the study area, with hourly model output. Both atmospheric and tidal forcing were included in the model, which used bathymetry sourced from Land Information New Zealand (https:// data. linz. govt. nz/ layer/ 51278-chart-nz-533-firth-of-thames/). Atmospher ic forcing was applied at a temporal resolution of 60 min and spatial resolution of 1 km. Due to the limited stratification in the FoT (Black et al. 2000, Zeldis et al. 2005), a 2-dimensional approach was used which was followed by a post-process transformation to obtain estimates of the depth-dependent currents. A logarithmic velocity profile was added to depth-averaged currents to simulate variations in horizontal water flow (Le Port et al. 2014). This high-resolution model was nested within a larger regional model of New Zealand-wide circulation at a grid size of 5 × 5 km. The high-resolution model was validated using tidal elevations obtained from a tide gauge at Tiritiri Matangi, located approximately 50 km northwest of the study site (36.605°S, 174.888°E) (Supplement 1).

Lagrangian particle tracking
Using the open source Lagrangian particle tracking software OpenDrift v.1.0 (Dagestad et al. 2018), we simulated individual dispersal trajectories for particles (representing recently spawned larvae) re leased from 5 selected aquaculture sites throughout the FoT (COR, MAN, WP, WB1, WB2; Fig. 1). Sites were selected to represent areas at which mussel aquaculture occurs in the FoT. Two sites were within (COR) or close to the entrance (MAN) of sheltered harbours; 2 sites represented the large exposed blocks of mussel farms in the eastern FoT (WB1 and WB2), and one site represented the small mussel farms in the western FoT (WP).  Okubo & Ebbesmeyer (1976) to account for subgrid-scale diffusion. A vertical diffusion constant of 0.01 m 2 s −1 over a 15 min timestep was included to account for turbulence in the fluid environment (based on parameters previously used for bivalve larval dispersal in Lundquist et al. 2004). A total of 4.9 million particles were released between the months of December and June, when peak settlement occurred over the years 1995−2005 (Alfaro et al. 2001, Smith 2019. Release sites were selected to represent all general areas in which mussel aquaculture occurs. Particles were evenly released within polygons delineated by mussel farming operations ( Fig. 1). As management and harvest practices are consistent within New Zealand (Jeffs et al. 1999) and therefore across the study area, the number of particles re leased from each site was scaled to the surface area occupied by each farm. Particles were released at a density of 10 particles km −2 of surface area of mussel farm at 2 h intervals to ensure particles were seeded across the tidal cycle. Full model parameters are presented in Table 1.
Particles were seeded at random depths within the water column between 1 and 10 m to simulate the depth of the dropper lines on which mussels are cultivated (Jeffs et al. 1999). No robust information exists on whether these mussels exhibit active horizontal or vertical swimming behaviour, though New Zealand bivalve larvae have been observed to be well distributed throughout the water column in a large harbour (Lundquist & Broekhuizen 2012). Particles were therefore given a terminal sinking velocity of 0.0025 m s −1 , which is at the lower range of settling velocities of larger bivalve larvae and may take into account some active settlement behaviour (Chia et al. 1984, Lundquist et al. 2009). Particles were released every 2 h for the full 7 mo settling period (Dec−June) every year from 1995−2005 and tracked for up to 5 wk. Age-dependent settlement was included in the model; if a particle encountered the seafloor or the coast (settlement habitat) between 3 and 5 wk after release (the settlement competency window) (Jeffs et al. 1999), it was deemed to have settled and remained at this location. Particles which did not en counter settlement habitat within this period were re tired from the model as this is beyond the known pelagic larval duration of this species (Hayden 1995). Due to the unavailability of data regarding the survival of mussel larvae, we assumed that mortality was spatially and temporally uniform. Thus, we did not include potential variations in daily mortality rates, but rather report proportional larval transport from each site.

Biophysical model analysis
Analysis of the OpenDrift output 'netcdf' files was conducted in MATLAB R2018a (MathWorks). Kernel density estimations (KDEs) were calculated using the methods of Botev et al. (2010). KDEs represented the probability that a released particle would settle in a given location. A total KDE was calculated for all particles released from all sites in the model. As the total KDE was dominated by particles released from areas with high surface area (WB1 and WB2), individual KDEs for particles released from each site were also calculated. The total proportion of particles which encountered suitable settlement habitat relative to the total number of particles released from each site over the 11 yr modelled period was also calculated.

Study sites
A total of 8 sites within the FoT were selected for trace elemental fingerprinting (C1−C7, W1; Fig. 1). Four of these sites were selected as they were within areas used for longline aquaculture and therefore potential larval sources. These included 3 of the modelled release sites (C3, C5, and W1). One of these aquaculture locations was situated in the western FoT at Waimang Point (W1), 2 were situated within Coromandel Harbour (C5 and C6), and 1 was situated within the large mussel farm blocks at Wilson Bay (C3). The remaining sites were selected as preliminary hydrodynamic modelling indicated they were within distances over which larvae from these aquaculture sites were able to disperse. Additionally, due to significant and varied anthropogenic modi fication of the study area over the last 100 yr (Hauraki Gulf Forum 2017), shell material produced at these sites is likely to differ in its trace elemental composition. Sites were generally between 2 and 5 km apart, although the largest distance between sites was 14.6 km (W1−C3) and the shortest was 1.5 km (C5−C6).

Field sampling
To determine natural variation in the trace elemental composition of mussels and establish reference trace elemental fingerprints within the study area, an in situ culturing technique was used (e.g. Becker et al. 2007, Kroll et al. 2018. Whilst previous studies have employed both moored and drifting in situ culturing methods to develop reference trace elemental fingerprints with differing results, it is difficult to establish which is more accurate (Kroll et al. 2018). We therefore elected to use a moored in situ culturing method for the development of reference trace elemental fingerprints. Macroalgae with juvenile mussels attached were harvested from the low intertidal zone at M ori Bay on the west coast of northern New Zealand (36.837° S, 174.427° E). This site was selected due to consistent year-round supply of juveniles and spat. Mussels were transported to the School of Biological Sciences at The University of Auckland, where individuals between 2 and 8 mm shell height were removed from the macroalgae and placed into in situ culturing containers (ICCs). Although incorporation of elements into shell material may vary with age (Strasser et al. 2008), previous research has shown that these ontogenetic differences in trace elemental fingerprints are unlikely to mask spatial differences (Norrie et al. 2016). ICCs were constructed from 750 ml food-grade, highdensity polyethylene jars with the lids and bases removed. Each end of the jar was covered with 500 µm Nitex ® mesh. This mesh size allowed for flow of water carrying nutrients but prevented the escape of animals. At least 20 juveniles were placed into each ICC, however up to 30 were translocated if enough animals were collected.
At each site, artificial settlement substrates (SSs) were deployed in order to collect recent settlers at each site. Plastic mesh tuffy scrubbing pads (SOS Tuffy, Clorox) were used as they have been extensively employed in previous settlement studies (reviewed by South 2016). Settlement on these substrates is widely accepted to be a proxy for larval supply (South 2016). These SSs were deployed with ICCs to form an experimental unit. To create an experimental unit, one ICC was attached to a rope at a depth of 7 m with one SS approximately 10 cm above the ICC and one 10 cm below. At each site, 3 of these experimental units were deployed on each sampling event (n SS = 3, n ICC = 3 site −1 per sampling event). Experimental units were deployed for ap proximately 5 wk, although inclement weather events resulted in different deployment lengths on each sampling occasion (ranging from 29−55 d) and the loss of all equipment at some sites for some dates. All sampling occurred over the austral summer and autumn (Jan−May) in a single year (2018), as this is a period of high settlement in the study area (Smith 2019). After retrieval, animals were removed from each ICC, placed into plastic zip lock bags, and frozen at −20°C until analysis. Each SS was placed into a zip lock bag in its entirety, as mussels are easier to remove from settlement substrate after freezing and thawing (Smith 2019).

Sample preparation for trace elemental analysis
Mussels from each ICC were defrosted and the shells were split open. One valve was randomly selected from each individual and any adhering particles were removed using stainless steel forceps. We inspected shells prior to translocation to ensure that shell material was a dark green. After recovery of ICCs, shell material with a distinct lighter green colour band following dark colourations was deemed to have been deposited in situ (Fig. 2a). The most recently formed portion of shell with a light green colour band along the axis of maximum growth was broken off and mounted onto microscope slides using double sided adhesive tape. If no distinct colour changes were observed, shells were not used for elemental analysis.
Settlers collected from 3 SSs at each site on each sampling occasion were analysed. In the laboratory, each SS was defrosted and all material was washed with tap water through a 200 µm sieve. All material washed off each SS was placed into deionised water (maintained at a resistivity of <1 MΩ cm), and mussels were sorted from this material. Mussels from each SS were selected and fixed to a microscope slide with double sided adhesive tape.

Trace element analysis
To determine the trace elemental composition of shell material, laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) was used. Elemental analyses were performed using a New Wave deep ultraviolet (193 nm) laser ablation system (Elemental Scientific Industries) coupled to an Agilent 7700 ICP-MS (Agilent Technologies). The laser operated with a spot size of 50 µm, a repetition rate of 5 Hz, and a dwell time of 30 s. Laser power was 45% and fluence was between 7 and 7.5 J cm −2 . National Institute of Standards and Tech nology (NIST) 610 and 612 glass standards were analysed every 20 spots for standardisation, calculation of internal precision, and calibration purposes. For full laser operating parameters and detection limits, see Supplement 2 (Tables S2 & S3).
All analyses were performed at the University of Auckland Plasma Mass Spectrometry Centre. Initially, 13 elements were monitored (Mn, Li, Co, Mg, B, Ca, Sr, Zn, Cu, Ti, Ni, Ba, Al). However, preliminary experiments revealed potential contamination from ICCs of Li, B, and Cu in shell cultured in these containers (Supplement 3), therefore these elements were excluded from further analyses. On reference shell material deposited in ICCs, one LA-ICP-MS spot was performed on the section of broken off shell approximately 200 µm from the ventral margin of shells (Fig. 2a). This material represented the most recently formed shell. On naturally settled individuals collected on the SS, the prodissoconch was analysed (Fig. 2b), as this reflects shell material formed at the individual's natal location (Becker et al. 2007, Carson et al. 2011. Data was processed using IOLITE trace elemental reduction software (Paton et al. 2011). Backgrounds were monitored for 30 s prior to each analysis. Data was then background-corrected by subtracting background average counts from the ablation counts. A pre-ablation procedure where the first 5 s of the ICP-MS signal was not included in the data was used to ensure that possible surface contamination or the periostracum was not included in the data (Marr et al. 2011, Norrie et al. 2016. To minimise the possibility of laser burn through to lower layers in the shell (Strasser et al. 2007), only the next 10 s of the laser dwell time was included in analyses. Data was then standardised using the most recent published NIST610 and NIST612 values (Jochum et al. 2005). All data was standardised to a trace element:calcium ratio (TE:Ca) in µmol of each TE to mol of Ca.

Statistical analyses
A discriminant function analysis (DFA) was performed on the TE:Ca ratios of shell material de posited in the ICCs. This allowed the examination of spatial variation in trace elemental fingerprints and development of site-specific reference signals. Due to differing covariances, a quadratic DFA (QDFA) was performed.  Elements were included stepwise into this QDFA to determine the optimal suite of elements for classification of shell material to its formation location (Dunphy et al. 2015, Norrie et al. 2016. QDFA statistical analyses were performed in JMP v.13 (SAS Institute). Due to the loss of sampling equipment at some sites through weather events over the study period, a QDFA was performed in which all samples collected from each site were binned together, regardless of the date collected. This allowed the development of mean trace elemental fingerprints for each site.
To predict natal origins, the discriminant function trained with the TE:Ca ratios of shell material de posited in ICCs was then used to assign the shell material at the umbo of each shell to a predicted formation location. As the loss of sampling equipment resulted in differing deployment durations for SSs, the number of settlers predicted to have originated from each site was standardised to a per day settlement rate. Receiver operating characteristic (ROC) curves were generated for each site to estimate the diagnostic ability of the QDFA. These ROC curves compared the true positive rate against the false positive rate. Posterior probabilities of each individual being correctly assigned to its formation location relative to other sites were also calculated in order to provide estimates of the confidence of the estimates.

Lagrangian KDEs
The total KDE of all particles released from all sites over the 11 yr modelled period (Fig. 3) indicates the area with the highest probability of receiving settlers was located on the eastern FoT to the north and south of the WB1 and WB2 release sites. Particles are likely to settle on the east and west coasts of the FoT (but not the south) with a higher likelihood of coastal settlement on the eastern FoT than the western FoT. The KDE indi-cated that benthic settlement is possible throughout the FoT, although it is concentrated in the east. The model also suggested a low probability of larvae being transported from the FoT to the greater Hauraki Gulf, particularly the area north of Coromandel Harbour and between Ponui Island and the mainland. This total KDE is dominated by particles released at areas with a high surface area from which more particles were released, therefore individual KDEs were generated for particles released from each site (Fig. 4).
The site-specific KDEs (Fig. 4) indicated that particles released from sites which were more sheltered Fig. 3. Kernel density estimates of all settled particles from all Perna canaliculus model release sites over all release dates in the Firth of Thames. Colour scale: probability that a particle will settle at a given location relative to the total number of parties released. Arrows: location of particle release locations that are difficult to see due to their size settled over a smaller area than those released from more exposed sites. This can be seen by particles released from COR exhibiting the highest probability of settlement in the relatively small area close to the mouth of Coromandel Harbour (Fig. 4a). Particles released from MAN also exhibited a high likelihood of settlement over a relatively small area close to their release site (Fig. 4b)  The 3 more exposed sites (WP, WB1, and WB2) had predicted settlement over wide areas throughout the FoT. Particles from WP (Fig. 4c) exhibited high settlement probabilities to the south and south east of the release site, including on the western coast of the FoT. A number of particles also settled on the eastern side of the FoT, particularly in the area south of Manaia Harbour. Particles released from the WB1 and WB2 sites (Fig. 4d,e) also settled throughout the FoT. The settlement patterns of particles released from these sites was broadly similar, with areas of high settlement probability located close to their release site as well as along the eastern coast of the FoT as far south as the town of Thames. Differences in in the settlement probabilities between the 2 sites were most evident in the south and northwest of the FoT. Particles released from WB1 were likely to settle in the southern FoT close to the Waihou River (Fig. 4d), whilst particles released from WB2 had a higher probability of settling in the north western FoT (Fig. 4e).

Spatial differences in the proportion of settled particles
Only 30% of the 4.9 million particles released between December and April across the years 1995−2005 encountered suitable settlement habitat (the coast or seafloor) during the settlement competency window. The remaining 70% of particles which did not settle were retired from the model. Differential settlement rates were observed between sites (Fig. 5). Highest settlement rates were predicted for particles released from the relatively sheltered sites COR (77% of particles settled) and MAN (53% settled). Conversely, the more exposed sites had a much lower proportion of settled particles which decreased as their distance from the coast increased: 35% of particles released from WP settled, 33% of those released from WB1 settled, and only 23% of particles released from WB2 encountered suitable settlement habitat during the settlement competency period.

Spatial variation in trace elemental fingerprints
Overall, the QDFA on reference shells correctly classified 72% of individuals to their growth location using 7 TE:Ca ratios (Mn:Ca, Co:Ca, Ba:Ca, Ti:Ca, Mg:Ca, Sr:Ca, Zn:Ca) ( Table 2). This is highlighted by the plot of the canonical scores from the QDFA (Fig. 6) indicating differences between elemental fingerprints of shell grown at different sites. However, individual site-level classification success was variable. High levels of classification success were observed at C2 (89%), C3 (100%), C5 (90%), and W1 (92%), indicating unique trace elemental fingerprints in shells deposited at these sites. Classification rates of shell material deposited at C1 was moderate, with 71% of individuals correctly classified. Approximately half of individuals grown at C4 (53%), C6 (55%), and C7 (52%) were correctly classified. The variability observed in the trace elemental fingerprints of shell material at C4, C6, and C7 shows that these sites are more difficult to differentiate based on their trace elemental fingerprint. The ROC curves were high for all sites (0.89−0.99), indicating the sensitivity of the QDFA was high (Table 2). Despite the relatively short distance between C5 and C6 (1.2 km), there were large differences in the trace elemental fingerprint of shell material deposited at these sites. Overall, 90% of individuals from C5 were correctly classified to their growth site. Classification success at C6, however, was only 55% with the majority of misclassified individuals attributed to C1, which was located furthest from C6 at approximately 16 km away. Release site Percent of released particles Fate of particle Retired Settled Fig. 5. Percent of Perna canaliculus particles released from each site which encounter suitable settlement habitat (defined as the seafloor or coast) during the settlement competency window 3−5 wk after release. Retired: the particle did not encounter suitable settlement substrate within the settlement competency window and was therefore removed from the model

Natal locations of recent settlers
The predicted larval origins of recent settlers collected at each of the sites monitored (Fig. 7, Table 3) suggested that the larval pool in the FoT is well mixed. Settlers collected at each site were a combination of larvae produced throughout the study region. A minimum of 5 and a maximum of 7 potential natal locations were predicted for settlers at each site. No settlers at any of the sites collected over the study period of January−May 2018 were predicted to have originated at C3. The confidence of the model in assigning individuals to their natal location was high, as shown by the mean posterior probabilities of correct assignment to each group ranging from 0.73−0.83 (Table 3) in agreement with the overall classification success rates of reference shell material (Table 2).
Although overall the larval pool was well mixed, some sites acted as larval sources and sinks relative to others. The most important larval source from the perspective of the highest number of larvae produced was C7. Of the individuals which settled at monitored sites, 30% were predicted to have originated from this location (1.63 settlers d −1 ) ( Table 3). The high contributions made by C7 were primarily driven by the high number of settled larvae originating from C7 at W1, which had the highest settlement rates of all sites monitored with 2.66 mussels d −1 settling at this location (Table 3). A large proportion of larvae which settled at C1 (48%; 0.28 settlers d −1 ) and C2 (46%; 0.25 settlers d −1 ) were also predicted to have originated from this location. A large proportion of larvae which settled throughout the study region were also predicted by the QDFA to have originated at W1 in the western FoT (0.95 settlers d −1 ) ( Table 3). Self-recruitment was generally low throughout the FoT with C7 (20% self-recruitment), W1 (14%), and C2 (15%) exhibiting the highest levels of the sites monitored. No settlers throughout the FoT were predicted to have originated at C3 despite the moderate settlement at this location. Only small numbers of larvae were predicted to have originated at both C2 and C5; although settlement at C5 was low, a moderate number of larvae settled at C2. These 3 sites can be considered larval sinks as few larvae from throughout the FoT were predicted to have originated at these locations.

DISCUSSION
This study set out to answer the following questions. (1) What is the potential for dispersal of Perna canaliculus from aquaculture locations throughout the FoT? (2) How do larval dispersal patterns generated by biophysical modelling and trace elemental fingerprinting differ? (3) What are the overall implications of any larval spill-over for restoration and aquaculture? The results of both the trace elemental fingerprinting study and the biophysical modelling indicated that there is potential for larval spill-over from aquaculture to areas throughout the FoT. The larval pool in the FoT is likely to be well mixed, with larvae from all sites potentially contributing settlers to, and receiving settlers from, other locations. Al though the results of both approaches showed similar broadscale results, different dispersal pathways were pre-dicted by each technique. These differences high light the importance of incorporating multiple tracking techniques for a holistic un der standing of larval dispersal patterns. Together, our results provide the first evidence, to our knowledge, that larvae produced in aquaculture can assist in the restoration of degraded wild bivalve populations through a larval subsidy).

What is the potential for dispersal of P. canaliculus from aquaculture locations throughout the FoT?
Both biophysical modelling and trace elemental fingerprinting indicated the Waimang Point mussel farms (WP, W1) were likely to be a source of mussel larvae which settle throughout the study area. In addition to contributing larvae to several sites, high settlement rates at this location also indicate that Waimang Point receives settlers from several locations. In a restoration context, the fact that this location may both receive and contribute larvae highlights the potential importance of this site in regional population dynamics. Restoration of benthic mussel beds at this site may provide habitat for larvae produced in other areas of the FoT as well as contribute larvae to other populations in the area.
Trace elemental fingerprinting predicted that the area to the north of Coromandel Harbour (C7) is an important larval source for settlers throughout the study area despite no known populations existing at this location. The biophysical model predicted settlement at this location was likely to be high, particularly for larvae released from Manaia and Coromandel Harbours (MAN and COR). One such explanation for this result may be that early stage larvae dispersing past this location incorporate elemental fingerprints from this site, although this is unlikely as hydrodynamic modelling indicated larvae in this area are retained close to their natal location for several days. Alternatively, these results present the possibility that a population exists at this location supported by aquaculture which may then in turn act as a source of larvae  (Table 3) for other sites in the FoT. This demonstrates the power of integrating the results of the 2 methods. The possibility of this stepping stone dispersal (Kimura & Weisss 1964, Palumbi 2003 highlights the need to consider system-wide connectivity into management actions based on source−sink dynamics (Kininmonth et al. 2011). Future surveys of potential mussel populations in the FoT should investigate this location to determine the presence of a population at this location which may supply larvae throughout the study area.
There are several larval sinks in the FoT study area. The area with the highest probability of settlement as predicted by the biophysical model was located to the south of Manaia Harbour, with larvae from each release site capable of settling at this location. Trace elemental fingerprinting also demonstrated that larvae settling in this area (C1, C2) likely originated from multiple sources. This finding is important as the area is known to support one of the few extant populations of mussels in the FoT (A. Jeffs, University of Auckland, pers. comm.), suggesting that aquaculture may already be supporting this population. The biophysical model also predicted that the area to the south of the Wilson Bay mussel farms would have high rates of settlement. Although this location currently supports no known populations of mussels, it is anecdotally known by local mussel farmers to be an area of high spat settlement (Smith 2019). Moderate settlement was also observed on artificial SSs in this area (C1, C2).
The similarity between the 2 methods of estimating larval dispersal shows that decisions on management in the study area based on this data can be made with a high degree of confidence and highlights the advantages of employing multiple techniques. The restoration of settlement habitat at these sink locations may allow the recovery of benthic mussel populations in this area due to high larval supply. Although adults are able to survive in this area (McLeod et al. 2012), it is important to consider factors which may prevent the survival of larvae or juveniles both during and after dispersal (Brandt et al. 2008, Pineda et al. 2010, Shima et al. 2015.

How do the 2 methods of estimating larval dispersal differ?
Although the overall pattern of a well-mixed larval pool in the FoT was demonstrated by both methods of tracking larvae, several differences in the pathways of dispersal were observed. These differences highlight the need to integrate multiple techniques of estimating larval dispersal and population connectivity to obtain a clearer picture of dispersal patterns (Ashford et al. 2010, Nolasco et al. 2018. The primary difference was the lack of settlers attributed to the large mussel farm blocks off Wilson Bay (WB1, WB2, C3) by trace elemental fingerprinting despite the prediction of the biophysical model that settlers  Table 3. Predicted formation locations of Perna canaliculus shell material at the prodissoconch of individuals which settled at monitored sites in the Firth of Thames between January and May 2018. As weather events resulted in non-standard deployment times for settlement substrates, standardised settlement per day was calculated. Settler contribution rate was calculated as the number of settlers per day of sampling which were predicted to have originated at each site from this location would dominate settlement patterns throughout the study area. The biophysical model also predicted that larvae produced within Coromandel Harbour (COR, C5, C6) would be retained in this area, although this was not shown to be the case in by trace elemental fingerprinting techniques. Several possible explanations exist for the differences between the 2 methods of tracking the dispersal of larvae in the FoT. It is likely that the variation in temporal scales over which the 2 studies were conducted was responsible for the differences observed. While the biophysical model considered circulation patterns over a 10 yr period, trace elemental fingerprinting only examined settlement patterns over a single spawning season. It is unlikely that dispersal patterns within a single year will fit with the 10 yr mean, as inter-annual variation in potential dis persal distances exists . This suggests that an extended sampling period should be undertaken to determine if long-term trends ob served through different techniques agree. In addition to possible variations in circulation, environmental conditions in areas of predicted high settlement may not have been suitable for the survival and settlement of larvae. The study area experienced heat wave conditions over the period in which field sampling was conducted (Jan−May 2018) (Herring et al. 2019), with reports of high mortality of mussels in aquaculture operations. This may have reduced both larval production at some locations as well as reduced the survival of larvae.
Additionally, there was variation in settlement rates predicted by the biophysical model between sites. These differences were driven by the rate at which larvae encountered potential settlement habitat within the time period in which settlement is possible. These results show that sites which produce the highest numbers of larvae may not necessarily contribute the highest numbers of settlers to populations. These differences in the rates at which larvae encounter settlement habitat have the potential to interact with direct stresses placed on an individual as it disperses, such as predation, as well as indirect dispersal stresses which may reduce post-settlement survival (Nanninga & Berumen 2014). Differences between the results of biophysical modelling and trace elemental fingerprinting highlight the need to consider realised dispersal and include empirical data into predictions of larval dispersal (Pineda et al. 2007). The similarities in overall patterns between the 2 methods employed in the current study highlight the utility of using computer-based simulations to guide the application of empirical validation experiments. These empirical experiments may then provide feedback towards validation and parameterisation of future models.

What are the implications overall for restoration and aquaculture?
The overall finding of larval spill-over from aquaculture has important consequences for the restoration and recovery of depleted bivalve populations. To our knowledge, this is the first study which demonstrates that bivalve aquaculture in coastal waters has the potential to provide a larval subsidy to relict and restored bivalve reefs and therefore assist with their recovery. This research shows that if bivalve aquaculture is conducted in an area with degraded bivalve populations, spill-over from these locations should be explicitly considered in restoration programmes. This research also shows that a network approach should be taken to restoration which includes connectivity between sites.
If larval spill-over from aquaculture to natural populations results in the recovery of degraded populations or assists with restoration efforts, connectivity to these locations should be maximised. This highlights the potential benefits of a network approach to restoration. The need for network approaches in marine reserve design is well known (e.g. Almany et al. 2007, Gaines et al. 2010, Green et al. 2015, Coleman et al. 2017, and these approaches should be translated to the design of restoration programmes. Where larval spill-over is desirable, additional bivalve aquaculture locations should be carefully planned in the context of the system in which this culturing takes place to ensure that larval spill-over to natural populations is enhanced. New restoration efforts should also be considered in the context of new or existing bivalve aquaculture to determine locations at which effort should be applied in order to maximise the chances of success (Arnold et al. 2017). The inclusion of aquaculture populations may simplify these network approaches, as many aquaculture populations are restocked after harvest resulting in relatively stable populations without the need to receive settlers.
It is also important to consider the possibility that larval export from aquaculture may negatively affect local populations of bivalves. The escape of finfish from farming operations has been shown to have potentially deleterious consequences on wild populations through the reduction of genetic diversity (Jørstad et al. 2008, Morris et al. 2008, Jensen et al. 2010). In bivalves, genetic introgression to wild stocks has also been demonstrated which may result in reduced fitness in the wild population (Apte et al. 2003). In cases where natural populations of a bivalve species exist, strategies which minimise potential larval export from aquaculture should be implemented. Potential strategies may include careful placement of aquaculture where settlement of larvae is unlikely, or the implementation of hatchery techniques such as spawning triploid bivalves with lower reproductive fitness (e.g. Honkoop 2003, Jouaux et al. 2010). In the case of the FoT, however, no large natural populations exist, making negative consequences less likely at present.

Future directions
This study has provided important insight into larval dispersal patterns from aquaculture and their interactions with restoration. However, to understand the ecological consequences of this dispersal it is essential that the survival to reproduction is quantified (Pineda et al. 2007, Nanninga & Berumen 2014. Several stressors may prevent the establishment of mussel reefs despite high larval supply (Burgess et al. 2012). In addition to spatio-temporal variation in mortality (White et al. 2014, Pineda & Reyns 2018, physiological stresses placed on an individual during dispersal, or post-settlement stresses may result in decreased fitness or lower reproductive output (Baker & Rao 2004, Burgess & Marshall 2011, Shima et al. 2015. Quantifying these stressors is essential for the selection and further assessment of viability of restoration sites, as larval supply and settlement will not result in a sustainable population if settlers do not survive. It is also important to quantify larval production in order to more accurately parameterise biophysical models and reconcile differences between modelled and empirically determined dispersal. For example, to better parameterise models, sites likely to receive high settlement as determined by modelling and empirical methods should now be monitored over long time scales in order to better understand inter-and intra-annual variation in larval production over the spatial extent of this study. Additionally, some locations are likely to result in lower larval output due to differences in the size structure of the population and age-dependent fecundity (Ren & Ross 2005). This is particularly important in an aquaculture setting where harvests occur at a relatively small size, which may reduce larval production relative to wild populations.
Understanding the consequences of varying dispersal distances is also important in the context of climate change. Changing ocean conditions such as salinity and temperature may result in changes in hydrodynamic flows on a number of scales (Wu et al. 2012, Sen Gupta et al. 2015. As these changes in hydrodynamic flows will directly affect the potential dispersal of larvae, it is important that changes in realised connectivity are considered (Coleman et al. 2011(Coleman et al. , 2017. Modelling exercises, for example, have shown that while oceanic circulation changes may affect larval transport, faster development and therefore earlier settlement (related to temperature) could override the effects of these changed transport pathways (Cetina-Heredia et al. 2015) or vice versa. It is also important to consider that faster development of larvae in warmer waters will only occur if the supply of nutrients is high enough to support this development; however, the distribution of plankton is also likely to change, which may decrease larval growth rates and even possibly increase mortality rates (Munday et al. 2009, Andrello et al. 2015. Finally, the possibility of employing methods such as those used by Nolasco et al. (2018), which explicitly incorporate the results of biophysical models and trace elemental fingerprinting into a single statistical framework, could prove very powerful. Acknowledgements. C.N. was supported by a UoA doctoral scholarship. Funding for fieldwork and laser analyses were provided by NIWA Coasts and Oceans core funding (COME1903). The hydrodynamic model was developed by MetOcean Solutions. M.R. was partially supported by the Moana Project (www.moanaproject.org), funded by the New Zealand Ministry of Business Innovation and Employment, contract number METO1801. Stuart Morrow from the UoA mass spectrometry centre aided with elemental analyses. Alex Vincent aided with figures. Jessica Bailey assisted in the laboratory. Fieldwork assistance was provided by Scott Edhouse, David Bremner, Peter Schlegel, Esther Stuck, Peter Browne, Errol Murray, and Jenny Hillman. Alan Bartrom of Gulf Mussels provided logistical support. Thanks to the Miller lab at OSU for their input on a draft of the manuscript. We also thank 2 anonymous reviewers and the editor for their constructive feedback which greatly improved previous versions of the manuscript.