Modelling long-term recruitment patterns of blue mussels Mytilus galloprovincialis : a biofouling pest of green-lipped mussel aquaculture in New Zealand

The green-lipped mussel Perna canaliculus forms the cornerstone of the New Zealand aquaculture industry. Like shellfish farming globally, P. canaliculus aquaculture is susceptible to the detrimental effects of biofouling. One of the greatest ongoing threats results from the recruitment of blue mussels Mytilus galloprovincialis onto cultured P. canaliculus and ropes used for collecting or on-growing juvenile stock. Using a dataset spanning 40 yr, we modelled spatio-temporal patterns in M. galloprovincialis recruitment in relation to a suite of potential explanatory variables and tested the ability of the best model to forecast recruitment. Our goal was to identify locations (sites and water depths) and times of peak M. galloprovincialis abundances, to enable the industry to implement management practices to minimise operational risks. Despite large inter-annual and spatial variability in recruitment patterns, our analyses revealed an upward trend in the abundance of M. galloprovincialis over the last 2 decades. Generally, seasonal patterns in abundance showed a large recruitment peak in October and a smaller one in April. There was a strong negative effect of depth, and weak effects of sea surface temperature and Southern Oscillation Index (in the previous month) on M. galloprovincialis abundances. The best model fitted the data well (R2 = 0.72); however, it provided only moderate forecasting power (R2 = 0.16), highlighting the challenges in forecasting recruitment. These results have been in corporated into a web application that enables aquaculture companies to interactively investigate historic trends in M. galloprovincialis and forecast abundances in the month ahead.


INTRODUCTION
Biofouling, the accumulation of organisms on submerged surfaces, is a significant operational problem for marine shellfish aquaculture globally (Adams et al. 2011, Fitridge et al. 2012).Suspended shellfish culture systems are particularly susceptible to the negative impacts of biofouling, in part because they provide a refuge from benthic predation and an ideal habitat for the proliferation of a range of problematic taxa.Biofouling can impact aquaculture operations at most production stages, including interference with spat (juvenile) recruitment, reduced growth, crop losses, decreased marketability and increased operational costs.
Spat are either wild-caught or hatchery-reared (Pérez Camacho et al. 1995, Lauzon-Guay et al. 2005, Díaz et al. 2014) and are a relatively vulnerable lifestage.As spat supply is at the start of the production chain, biofouling of spat can have significant 'downstream' effects on industry profitability.Since the cost and logistic barriers to managing problematic biofouling species after they have established on aquaculture installations are considerable, there is increasing global interest in preventative strategies that seek to minimise colonisation by pest species, not only at the spat settlement/recruitment stage, but also at the line seeding and crop grow-out stages of industry operations (Fitridge et al. 2012, Fletcher et al. 2013, Sievers et al. 2014).
Implementation of preventative approaches, for example avoiding or minimising peaks in colonisation by problematic species during key stages of production, requires knowledge of spatio-temporal patterns of their occurrence, and of the factors that determine their recruitment.This information could then be used to facilitate the development of prevention-based management measures.The timing and magnitude of recruitment may be controlled by various factors, including larval supply (Cáceres-Martínez & Figueras 1998a,b, Hoffmann et al. 2012), the presence of conspecifics (McGrath et al. 1988), water depth (Fuentes & Molares 1994) and environmental variables such as water temperature (Avendano & Cantillanez 2013), food availability (Toupoint et al. 2012), wind patterns (Barria et al. 2012), waveexposure and water flow (Hunt & Scheibling 1996), and solar irradiance (Fuentes-Santos et al. 2016).For example, Fletcher et al. (2013) considered seasonal recruitment patterns of the problematic ascidian Didemnum vexillum in relation to sea surface temperature (SST), to identify low-risk 'windows' for undertaking critical mussel industry activities such as spat deployment to nursery areas.
In New Zealand, biofouling is emerging as one of the significant constraints to the development of the mussel-farming industry, which is based on cultivation of the indigenous green-lipped mussel Perna canaliculus (hereafter Perna).In the country's most significant aquaculture region at the north end of the South Island (Fig. 1A), the native New Zealand blue mussel Mytilus galloprovincialis (hereafter My tilus) has been a persistent and wide-spread biofouling problem for several de cades.This species is part of the Mytilus edulis complex (Westfall & Gardner 2010, Oyarzún et al. 2016).It represents a Southern Hemisphere lineage of the Mediterranean mussel Mytilus galloprovincialis, with evidence of extensive hybridisation and backcrossing between New Zea land blue mussels and the invasive Northern Hemisphere lineage of M. galloprovincialis (Gardner et al. 2016).The Northern Hemi sphere lineage is considered an invader in many parts of the world with the ability to rapidly become established (Apte et al. 2000), although it is also an important aquaculture species in many countries (FAO 2004(FAO −2016)).In New Zealand, Mytilus can pre-empt space on Perna spat-catching ropes, and recruit directly onto Perna spat (a process referred to by industry as 'over-settlement'), thereby significantly reducing spat supply.For example, 104 Mytilus that is transferred onto nursery ropes during Perna seeding processes can significantly decrease Perna retention on the lines, with Perna seed losses in some cases exceeding 95% (Carton et al. 2007).
Mytilus also has a range of additional impacts at subsequent stages of production, including causing significant losses during the final crop grow-out stage (Forrest & Atalah 2017).
In this paper, we used spatio-temporal Bayesian modelling to analyse a long-term regional dataset on patterns of Mytilus recruitment to settlement arrays deployed throughout the north of the South Island region, and relate patterns to a suite of potential explanatory variables.Specific goals were to: (1) identify times and locations of greatest Mytilus recruitment risk, (2) identify environmental drivers and other explanatory factors affecting recruitment patterns and (3) illustrate how such information can be used to develop mitigation measures that can be incorporated into husbandry practices.As part of the latter, we describe the development of an interactive data mining and forecasting tool that is of relevance to biofouling management in aquaculture internationally.

Mytilus spat data
Our study used datasets obtained from the New Zealand Marine Farming Association which contained mussel spat monitoring data dating back to 1975.Over this period of ca.40 yr, monitoring of Mytilus has been conducted alongside monitoring of Perna using spat collectors consisting of 0.5 m sections of mussel industry 'Xmas tree' rope (Fig. 1B).Monitoring is ongoing, and involves deploying the collectors for 1 wk periods at depths from 1 to 15 m, at multiple sites across the study area.These sites represent some key mussel aquaculture sub-regions, namely Pelorus Sound, Queen Charlotte Sound, Port Underwood, Tasman Bay and Golden Bay (Fig. 1A).After 1 wk, the spat collectors are retrieved and taken to a laboratory where Mytilus spat are counted using a binocular microscope.At this stage the recruits are ca.0.5 mm in shell length.

Environmental data
Site-specific SST satellite data were obtained from the Group for High Resolution Sea Surface Tempera-ture (https://www.ghrsst.org/).The data covered the time period from January 2002 until December 2015.Southern Oscillation Index (SOI, Trenberth 1984) data were ob tained from the Department of Science, Information Technology and Innovation, Queensland, Australia (https://www.longpaddock.qld.gov.au).The SOI is a measure of the large-scale pressure differential be tween the western and eastern tropical Pacific coinciding with El Niño and La Niña episodes (Troup 1965).In general, smoothed time-series of the SOI correspond well with changes in ocean temperatures across the Pacific.

Statistical analysis
Initial data exploration revealed geographic and temporal gaps in the dataset, especially in the earlier years of the monitoring programme.As such, while we describe trends in Mytilus since 1975, the statistical analysis focused on a subset of data covering a continuous 14 yr period from 2002 to 2015 (Table A1 in the Appendix).The dataset from that pe riod was relatively complete in terms of seasonal, geographical and water depth representation, and consisted of 19 968 observations of Mytilus spat count m -1 of rope (the response variable).The dataset also included a suite of potential continuous explanatory variables (year, month, depth, latitude/ longitude, SOI, SST), which were used to model Mytilus spat abundances in space and time.The model was constructed using a Bayesian approach based on integrated nested Laplace approximations (INLA, Rue et al. 2009).This methodology is freely available in the statistics software R (R Core Team 2014) as part of the package R-INLA (www.r-inla.org).This package has become popular among applied research in a wide range of fields, including ecology, criminology and epidemiology (Blangiardo & Cameletti 2015).Bay esian methods provide a power ful and flexible approach for large or complex data sets, such as the Mytilus spat dataset described here.
The model was built by approximating the predicted distribution of spat m -1 given previously observed Mytilus occurrences and their associated explanatory variables as described above.As the raw data consisted of a large number of zero counts for Mytilus (55%) and high over-dispersion, we aggregated the total Mytilus spat m −1 to Year-Month-Site-Depth nodes (Fig. 1C) and included an 'offset' term to account for the number of samples that could have contributed to this total.This approach reduced the space and time dimensions of the model and allowed for more accurate approximations of the true probability distribution.
Our base model specified the covariates SOI and SST, lagged by 1 mo to imitate forecasting conditions, as linear fixed effects and water depth as a non-linear random effect.Long-term effects and inter-annual variability were included as spatiotemporal random effects, where: (1) the spatial effects were linked to the spatial random field at each mesh node (Fig. 1C), and (2) each month's spatial effect was estimated under the assumption that the spatial pattern changed each month and was correlated with the spatial effects of the previous and subsequent month.The structure of this model allowed us to infer Mytilus spat abundances at previously sampled sites as well as unobserved locations.
The response variable, total monthly Mytilus spat count m −1 , was assumed to follow a negative binomial distribution, with the parameters of this distribution (φ) linked to a structured additive predictor η through a link function μ = E(e η ), where log(E) is equal to the offset of η.The linear predictor η for the final model was defined as: (1) where β 1 and β 2 are the linear regression coefficients for the fixed effects, lag 1 (i.e.previous month) SOI and SST, respectively.The third term is the sum of random walk order 1 smooth functions defining the random effect of depth where coefficients vary depending on the observed depth values, w k .The final term is the spatio-temporal random effect, with Matérn correlation structure de fining the spatial effect, and additional autoregressive (order 1) correlation structures for consecutive months.
The strength of the association between observed counts of Mytilus spat versus fitted counts from the historical dataset was assessed using Pearson correlation.To test the power of the final model to forecast Mytilus spat abundances 1 mo ahead in relation to current environmental conditions at a given location and time, we used 1 mo out-of-sample cross-validation.This technique uses environmental conditions for a given month (i.e.SOI and SST) to forecast Mytilus abundance 1 mo ahead.This process was repeated for the previous 24 mo to capture any seasonal variation.Model performance was checked by plotting ob served versus both fitted and forecasted values along with the 1:1 line.Additionally, model performance was summarised using 3 statistics; regression R 2 , Nash-Sutcliffe efficiency (NSE) and percent bias (PBIAS) for both fitted and forecasted values.The re gression R 2 value is the coefficient of determination derived from a regression of the observations against the predictions.NSE indicates how closely the observations coincide with predictions (Nash & Sutcliffe 1970).NSE values range from −∞ to 1.An NSE of 1 corresponds to a perfect match between predictions and the observations.An NSE of 0 indicates the model is only as accurate as the mean of the observed data, and values less than 0 indicate the model predictions are less accurate than using the mean of the observed data.Bias measures the average tendency of the predicted values to be larger or smaller than the observed values.Optimal bias is 0, whereas positive values indicate underestimation bias, and negative values indicate overestimation bias (Piñeiro et al. 2008).PBIAS is computed as the sum of the differences between the observations and predictions di vided by the sum of the observations (Moriasi et al. 2007).A rule of thumb is that model predictions are satisfactory if NSE > 0.50 and if PBIAS < ± 25%, and are good if 0.65 < NSE > 0.75 and if 25% < PBIAS < ± 40% (Moriasi et al. 2007).Additionally, we used ac cepted conventions to interpret correlation strength according to R 2 values (Cohen 1988, Evans 1996).

RESULTS
Over the 40 yr period represented by the dataset, there was a large inter-annual variability in the recruitment of Mytilus across the study area (Fig. 2A).Low abundances of spat (mean < 103 ind.m −1 ) were recorded until the early 1990s, with the exception of the first year of monitoring in 1975 (mean ± SE = 330 ± 105 ind.m −1 ).However, after 1994, there was a general trend for a steady annual increase in the abundance of Mytilus spat, with peaks in mean abundance of 464 ± 145 and 797 ± 197 m −1 recorded in 1995 and 2005, respectively.The upward trend in the last 2 decades is described by the smoother timeseries in Fig. 2A.
The final model consisted of the fixed effect SOI t -1 , SST t -1 , the non-linear effect of depth and the spatiotemporal effect (see 'Materials and methods').Overall, the recruitment of Mytilus throughout the year had 2 distinctive peaks: a larger one in the Austral spring through to early summer (i.e. between October and December) and a smaller one in early autumn (i.e.April, Fig. 2B).However, at a finer resolution, the model predicted distinct geographic patterns in the distribution and abundance of Mytilus, which varied considerably across months and seasons (Fig. 3).At sites in the western sub-region (i.e.Golden and Tasman Bay) the autumn peak extended until May and was of comparable or even greater magnitude to the spring peak (Fig. 3).At sites in the easternmost sub-region (i.e.Queen Charlotte Sound) a strong early autumn peak was not evident, with the model revealing greatest spat abundances in June and October (Fig. 3).Generally, there were more pronounced peaks in mid-Pelorus Sound, compared to the other sub-regions (Fig. 3).Increasing water depths were generally associated with lower expected Mytilus abundances (Fig. 4).There was a positive effect of water depth on spat abundances until 6 m, and a negative effect of depth between 6 and 15 m (Fig. 4).There was a marginal negative effect of SST (mean effect = −0.077,95% CI = −0.312 to 0.157), indicating that lower Mytilus abundances were correlated with higher SST in the previous month.There was a small positive effect of SOI on spat abundances (mean effect = 0.009, 95% CI = −0.024 to 0.042).Although neither SST nor SOI were deemed significant based on 95% confidence intervals, they were retained in the final model, as cross-validation indicated that they increased the model's explanatory and forecasting power.Overall, the validation process indicated that the final model performance was good (Fig. 5A, R 2 = 0.72 and NSE = 0.72), and there was a minimal overestimation of 0.9% as indicated by the PBIAS statistic.On the other hand, under forecast conditions (1 mo ahead), the model performance was only moderate (R 2 = 0.16 and NSE = −0.33,Fig. 5) and over-predicted Mytilus spat abundance considerably, particularly at low abundances (PBIAS = 82%, Fig. 5B).

DISCUSSION
Biofouling poses a significant threat to shellfish aquaculture globally and may drastically affect spat (juvenile) supply, crop yield and marketability, and industry operations (Adams et al. 2011, Fitridge et al. 2012).Shellfish aquaculture that relies on wildcaught spat is particularly vulnerable to the negative impacts of biofouling.Given the lack of cost-effective management options to mitigate the impacts of established biofouling, a potentially feasible approach is the implementation of husbandry practices to minimise exposure to initial colonisation by problematic species (Sievers et al. 2014).Here we described high variability in spatio-temporal patterns in a regional dataset on Mytilus recruitment spanning 40 yr.A Bayesian spatio-temporal modelling approach was applied to a 14 yr subset of the data, enabling identification of locations and windows of time when Perna farms have been historically prone to colonisation by Mytilus.The model was simultaneously evaluated for its ability to forecast Mytilus spat abundances 1 mo ahead.

Spatio-temporal patterns
Large inter-annual variability in abundances, such as evident in the data analysed here, is an intrinsic characteristic of the recruitment of mussels (Hunt & Scheibling 1998, Porri et al. 2006) and other marine invertebrates (Rodríguez et al. 1993).Processes that cause such variation may include fluctuations in physical factors, resource availability and biological interactions (Underwood & Chapman 2000).In this Although the cause of the apparent long-term increase in recruitment cannot be inferred from this study, a possible explanation relates to the considerable expansion of mussel cultivation since it began in the 1970s.Overall, there has been an increase in the mussel farming area, and hence habitat for Mytilus, from ca. 100 ha in the mid-to late-1970s to ca. 2500 ha present day (Handley 2015).In fact, in Pelorus Sound (where the majority of farms in our study region are located) there was a moderate positive correlation between the cumulative area of mussel farms and the annual median Mytilus abundance (r = 0.63, p < 0.001).An in creased Mytilus population could lead to an increased larval inoculation pressure for ongoing establishment (Ruiz et al. 2000, Lockwood et al. 2009), such that increased abundances on mussel farms could become a self-perpetuating problem.A similar phenomenon was suggested for salmon farms in Norway, where the recruitment of Ectopleura larynx and other sessile biofouling species was dramatically enhanced within farms relative to control locations (Bloecher et al. 2015).
Although Mytilus recruitment was recorded through out the year, there was a strong overall seasonal effect, with 2 peaks: a major and a minor peak in October and April, respectively.Seasonal fluctuations in Mytilus recruitment of comparable magnitude have been previously reported in New Zealand (Kennedy 1977, Meredyth-Young & Jenkins 1978) and other parts of the world (Dix & Ferguson 1984, Yildiz & Berber 2010, Yildiz et al. 2010).These fluctuations are most likely related to the reproductive seasonality observed in this species, with a peak spawning season in spring coinciding with increases in water temperature and chlorophyll a concentrations.Several spawning events can occur after the first one, generally until summer or late autumn (Cáceres-Martínez & Figueras 1998b).In autumn, a second peak in food availability (Gibbs et al. 1992, Gibbs & Vant 1997, Ogilvie et al. 2000, MacKenzie & Adamson 2004) may favour the accumulation of reserves for gametogenesis and a subsequent minor spawning event.
When the long-term and seasonal patterns in Mytilus abundance were considered by the main sub-regions, distinct geographic variation was evident (see Fig. 3).For example, Golden Bay sites showed a minor end of summer peak between April and May, rather than between March and April as observed in Pelorus Sound.The greatest Mytilus abundances overall occurred in mid-Pelorus Sound and Queen Charlotte Sound, during October and July.The mid-Pelorus Sound is the most intensely farmed area of the study region (Aquaculture New Zealand 2012).As discussed above, Mytilus populations established on these farms may act as a larval source and enhance subsequent recruitment in adjacent areas.Additionally, mid-Pelorus is a relatively enclosed area, with several inner bays having low water flow (< 0.1 m s −1 ) and slow flushing times (Gibbs et al. 1991), which may favour larval retention and settlement (McQuaid & Phillips 2000).
The overall pattern of higher Mytilus recruitment in shallower (< 6 m) compared to deeper waters (> 6 m) is consistent with the depth patterns reported in previous studies in Pelorus Sound (Meredyth-Young & Jenkins 1978) and elsewhere (Cáceres-Martínez & Figueras 1997, Holthuis et al. 2015).The occurrence of lower Mytilus spat abundances in deeper waters supports the existing practice of some mussel farmers, who submerge lines to avoid or minimise over-settlement by Mytilus.These depth differences are conceivably explained by the vertical distribution of Mytilus larvae in relation to environmental variables.Competent larvae of Mytilus are generally not homogeneously distributed across the water column and have a tendency to aggregate in surface waters above the thermocline/halocline (Dobretsov & Miron 2001).Mussel cultivation waters across our study region experience a thermocline at ca. 10 m or deeper, and in Pelorus Sound surface salinity at shallow depths <10 m can be reduced due to significant freshwater inputs (Gibbs et al. 1991, Gibbs 2001).

Utility and operational implications of the model
Relationships between Mytilus spat recruitment and environmental conditions, such as SST and SOI, were generally weak but did at least provide additional accuracy to the model.The marginal effect of SST (in the previous month) was not as great as might be anticipated given the likely importance of temperature in the reproductive seasonality of Mytilus noted above.More broadly, water temperature is generally recognized as a driving factor in the settlement and recruitment of marine invertebrates, as it can affect adult reproductive output, and the subsequent survival, growth and metamorphosis of larvae (Kinne et al. 1970, Coma et al. 2000, Bates 2005).
The significant spatial random component in the model likely reflects the effect of other unmeasured factors that influence the recruitment patterns observed.As well as the probable importance of chlorophyll a as suggested above, other environmental factors such as salinity, wind patterns and hydrographic characteristics may influence reproductive seasonality and recruitment.Additionally, a very recent study of settlement patterns in Mytilus highlighted solar irradiance as a key predictor variable, suggesting that solar irradiance during late winter indirectly drives the timing and intensity of settlement onset (Fuentes-Santos et al. 2016).Recognising that it can be difficult to attribute changes in biological responses to a single variable (Coma et al. 2000), it is clear that a broad range of variables would ideally be incorporated into our model.Predictive power may have increased by incorporating measures of mussel food sources, such as phytoplankton and detritus, or other environmental variables (e.g.salinity and wind).For example, the recruitment of many coastal species with long larval duration, such as Mytilus (2−4 wk) is driven by surface currents, which are largely controlled by wind patterns (Rough garden et al. 1988).Such advances may be possible in future years, especially with the development of monitoring technologies that facilitate the collection of large and long-term environmental datasets (He et al. 2015).Even with comprehensive environmental data, however, recruitment variability is likely to be attributable, in part, to unrelated factors, including attributes that are intrinsic to the mussels themselves.For example, mussels have the capacity for secondary settlement by attaching and detaching from the substrate several times until they find a suitable substrate (Bayne 1964, Newell et al. 2010).In the case of Mytilus, the process of secondary settlement may occur for several weeks after initial larval settlement (Bownes & McQuaid 2009).As such, the recruitment response patterns described in the present study, and indeed other studies of early post-larval mussel settlement and recruitment, may not reflect the final distribution patterns of adult mussels (Le Corre et al. 2013).
Despite the highly variable nature of the data, and the weak individual effects of potential environmental drivers, the model nonetheless provided relatively reliable estimates of Mytilus spat abundances (observed vs. fitted data, R 2 = 0.72).These findings can be used to guide farm management practices, such as the location and timing of deployment of Perna spat collection ropes to minimise Mytilus over-settlement.However, the moderate relationship between observed and forecasted data (R 2 = 0.16) suggests that the prediction of Mytilus recruitment 1 mo ahead should be viewed with caution.Clearly, there are phenomena contributing to recruitment patterns that have not been captured by the model and underlying data.We can nonetheless provide a reasonable model of the likely general abundance of Mytilus spat at given locations and times.
Considering that the regional economic loss of Mytilus on Perna aquaculture has been estimated at ca. US$16 million annually in the Marlborough Sounds region alone (Forrest & Atalah 2017), the ability to implement management practices to avoid high abundances of Mytilus may translate into considerable economic benefits.Predictive maps such as depicted in Fig. 3 provide a useful means to illustrate to marine farmers the risk periods and hotspots for Mytilus recruitment.For example, a Mytilus recruit count of > 200 m −1 , evident in orange and red colours in Fig. 3, is the approximate threshold that the industry considers operationally significant.This information can be considered in relation to peak times of Perna settlement in the region (i.e.April, Atalah et al. 2016).Accordingly, the spatio-temporal model has recently been integrated into a web application (https:// cawthron.shinyapps.io/BMOP/)described by Atalah et al. (2016).The web application currently displays spat count data for both Mytilus and Perna, for more than 50 monitoring sites across the study region.The application enables marine farmers to interactively explore historic settlement patterns of Mytilus and Perna spat at fine-scale resolution (i.e. at particular farm sites or bays), with a future step being to extend the forecasting model developed for Mytilus to also include Perna.This tool will therefore facilitate management decisions that aim to maximise Perna collection while minimising losses due to Mytilus (e.g. based on the ratio between Perna and Mytilus spat abundance), and will therefore provide an integrated production and pest management tool for famers.This type of application is clearly of relevance to the management of any fouling pest in aquaculture for which spatio-temporal variation in occurrence provides opportunities for management (Fletcher et al. 2013, Sievers et al. 2014).Although the tool does not provide a panacea for the management of all high-risk biofouling, by addressing the most significant species it facilitates sustainable and efficient aquaculture production.

Fig. 1 .
Fig. 1. (A) Study area, in relation to New Zealand's South Island (inset), showing the main sub-regions.(B) 'Xmas tree' rope used in collectors to sample Mytilus spat.(C) Mesh pattern used to calculate the effect of geographical location (Gaussian random field) and the location of the study sites (red dots)

Fig. 2 .
Fig. 2. (A) Long-term and (B) seasonal trend in mean Mytilus spat abundance m -1 rope.The red dashed line is the smoothed time series using a cubic spline and 95% confidence intervals (grey shading)