From individual vital rates to population dynamics: An integral projection model for European native oysters in a marine protected area

Funding information Natural England, Grant/Award Number: SR17/ 18‐1031455‐008; Natural Environment Research Council, Grant/Award Number: NE/L002582/1 Abstract 1. Following an 85% decline in global oyster populations, there has been a recent resurgence in interest in the restoration of the European native oyster Ostrea edulis. Motivations for restoration from environmental stakeholders most often include recovering lost habitats and associated biodiversity and supporting ecosystem function. In coastal communities, another important justification is recovery of traditional and low‐impact fisheries but this has received less attention. 2. Many restoration projects across Europe focus on the translocation of adult stocks, under the assumption that the limit to population growth and recovery is adult growth and survival. This may not necessarily be the case, especially where knowledge of large extant adult populations exists as in the Blackwater, Crouch, Roach and Colne Marine Conservation Zone in Essex, UK. Identifying what limits population growth for restoration and recovery is an important conservation tool. 3. Here, the first size‐dependent survival, growth and fecundity data for free‐living O. edulis from a novel field experiment are used to parameterize an Integral Projection Model that examines the sensitivity of a flat oyster population to variation in individual vital rates and to potential harvesting – an original objective of a coastal community‐led restoration project. 4. Given the high adult fecundity in this species, population recovery is most sensitive to changes in recruitment success; however, elasticity (proportional sensitivity of the population) is more evenly spread across other parameters when recruitment is already high. Based on locally agreed management objectives, recovery to double the current stock biomass should take 16–66 years (mean = 30 years) without active intervention. At that point, harvest rates could be sustained below 5% of the harvestable adult size whilst ensuring λs remains above 1.


| INTRODUCTION
Once a common species around all coastlines of the UK and widespread throughout Europe, the European native oyster, Ostrea edulis, has now been reduced to a few remaining strongholds, with oyster populations of all species thought to have declined by as much as 85% worldwide (Beck et al., 2011). In recent years, however, there has been a resurgence in European native oyster restoration across the historical range of this species (Pogoda, Brown, Hancock, & von Nordheim, 2017).
Many oyster restoration projects cite an ecosystem service as their primary motivation for restoration. These range from potential biodiversity benefits from the presence of flat oyster beds or reefs (Lown, 2019;Pogoda, 2019) to potential filtration benefits resulting in increased water quality (Rodriguez-Perez et al., 2019;Wilson, 1983).
While not yet common in oyster restoration in Europe there is also a long history of restoration and conservation interventions for restoring traditional shellfish fisheries (Allison, Hardy, Hayward, Cameron, & Underwood, 2020).
Many European native oyster restoration projects are directed towards the translocation of adult stock (i.e. the movement of oysters from coastal areas of high/higher abundance to areas of low abundance), with population status largely assessed via estimated abundance or biomass of adults (Marine Management Organisation, 2019). This is in part due to the very low abundance or extinction of native oysters from many areas (Fariñas-Franco et al., 2018). However, this may not always be the most appropriate method to assess and restore O. edulis populations, especially in those areas where oysters already have significant local stocks (e.g. Scotland, Ireland, Essex) (Eagling, Ashton, & Eagle, 2015;Lown, 2019;McGonigle, Jordan, & Scott, 2016). It is also good conservation practice to a priori examine what life history stages are most likely to respond to restoration or management interventions, and thus affect population growth (Montero-Serra et al., 2018). Such comprehensive assessments of population status and sensitivity require real-time measurements of demographic parameters such as size structure, growth, survival and fecundity of individuals in the population (Caswell, 1989;Ellner & Rees, 2006;Merow et al., 2014).
Once undertaken, these modelling approaches can determine if populations are likely to be influenced by changes to recruitment or survival, thus informing management. For example, increasing spat settlement and survival via manipulation of habitat may be a more appropriate restoration tool than translocating adult oysters, or halting their fishing in a given area.
Previous studies have monitored the growth rates and survival of various oyster species in aquaculture or laboratory settings in controlled environments, in bags at high density above the seabed or in cages. However, data in a useable format for modelling are rarely published (Katkansky, Dahlstrom, & Warner, 1969;Pogoda, Buck, & Hagen, 2011). In addition, no study to date has simultaneously monitored individual vital rates such as growth and survival of flat oysters at naturally occurring densities (Allison, 2017;Helmer et al., 2019;Lown, 2019), particularly on the sea bed where smothering, predation risk and food resources may differ from conditions when raised above the benthos or in aquaculture settings (Sawusdee, 2015;Zwerschke, Emmerson, Roberts, & O'Connor, 2016). In order to appropriately monitor the success of restoration projects it is essential for these demographic data to be regularly estimated in-situ and incorporated into model predictions of responses to restoration activity.
A number of different types of population model are available to assess and make predictions on populations using vital rates such as growth, survival and fecundity. Individual based models follow every individual within a population to assess individual outcomes and consider individual-specific characteristics and dynamics, with population dynamics a sum of these outcomes. Conversely, distribution based models follow populations and their dynamics via population-level distributional changes (Picard & Liang, 2014). Subsequent differences between model types are largely due to the continuous or discrete nature in which age or size, time and/or reproduction are treated. Examples of these range from matrix based models whereby both time and size or age are discretized, to Ordinary Differential Equations or Physiologically Structured Population Models where time, reproduction and size are usually treated as continuous variables (de Roos & Persson, 2013;Picard & Liang, 2014).
Integral projection models (IPM) also use discretized time, calculating the dynamics of, and changes in, abundance and vital rate distributions of populations over fixed periods of time. Time segments can be set to any unit, e.g. 1 day, 1 month or 1 year, depending on the lifespan of the species under investigation (Merow et al., 2014). IPMs use a series of regression equations to parameterize growth, survival and reproduction rates to incorporate individual-based variation within the growth transition part of the model (Picard & Liang, 2014).
IPMs use an integral equation called a kernel to describe changes of state of individuals from one timestep to another (e.g. their growth, survival and fecundity and how these are likely to change with each other throughout an individual's life; Merow et al., 2014;Rees, Childs, & Ellner, 2014). This means that IPMs are useful in situations where abundance estimates from census data are based on discrete timesteps whereas data on life stage-dependent vital rates have been collected as a continuous distribution, such as annual growth rates and survival. For this reason, an IPM of an O. edulis population and population dynamics is presented with how it can be used in restoration projects.
This study is of the native oyster populations in the Blackwater, Crouch, Roach and Colne Estuaries Marine Conservation Zone (BCRC MCZ) in Essex in the southern North Sea. The overall objective of the MCZ is to protect and recover native oysters and their habitats, with the population maintained in numbers which enable it to thrive (UK GOV, 2013). Undertaking active intervention is recognized as required to achieve this aim (Allison, 2017;Helmer et al., 2019). The Essex Estuaries have a long cultural history of oyster fishing, dating back to Roman times, with the Colne and Blackwater oyster fisheries listed in the Domesday Book (Benham, French, & Leather, 1993 The aim of this study is therefore to obtain locally relevant demographic information for native oysters in the BCRC MCZ and use this in an IPM modelling framework to assess how long recovery of the population may take under current growth, survival and reproduction rates and, subsequently, how a 'restored' population may respond to reintroduction of harvestingone of several restoration objectives.

| The life cycle of O. edulis
The European native oyster is a protandric, sequentially hermaphroditic, slow growing bivalve. It spawns in the summer months (June to September) with offspring experiencing up to 10 days of brooding then up to 10 days of larval phase before settling as spat onto hard substrate (Helmer et al., 2019). Unlike the more widely studied Crassostrea spp., O. edulis is largely subtidal and is ovoviviparous. In addition, O. edulis are notoriously difficult to age, not exhibiting clear nacreous rings on the exterior of the shell. When these rings are visible, O. edulis may occasionally lay down multiple rings in a single growing season (Orton & Amirthalingam, 1927;Richardson, Collis, Ekaratnc, Dare, & Key, 1993). This means that IPMs incorporating age, as have been created with Eastern oysters (Moore, Lipcius, Puckett, & Schreiber, 2016;Moore, Puckett, & Schreiber, 2018), are not currently appropriate for this species.

| The IPM
An IPM maps the distribution of some form of vital rate that changes over an individual's life, such as size or age at time t, to the distribution of that same vital rate at time t + 1 (Merow et al., 2014). For the model of O. edulis in Essex, the main vital rate used to predict changes in the population was body size, data that is collected by most fisheries and restoration programmes globally. Area of the oyster (mm 2 ) was used as a proxy for body size assuming each oyster to have an elliptical shape [i.e. area = (height/2)(length/2)π]. This was deemed most appropriate owing to the nature of oyster growth, with some oysters growing in length rather than shell height (Lown, 2019). This model was developed following the collection of three primary types of data. These data were growth and survival data from an individualbased monitoring experiment (hereafter referred to as the string experiment), census data from multi-year extensive dredge surveys, and finally, reproduction and fecundity data from a literature review and assessment of the prevalence of fecund females during the breeding season. All field data were collected from within the BCRC MCZ between 2014 and 2019, during which time no harvesting of oysters was permitted to occur. ing unknown values through interpolation within a specified area, weighted by the value of points closest to those being estimated (Chen et al., 2016). A power value of 2 was used with five points used to calculate each interpolated cell owing to the five-point sampling design of the census survey and to give higher power of influence of near sites over more distant ones. An output cell size of 10 × 10 m was used to speed processing whilst maintaining a highlevel detail within the calculation. Data from post-winter (March) surveys were used for the IPM census data.

| String experiment: individual-based data
To monitor growth and survival rates of individual oysters, an in situ experiment was designed. This was largely based around tile-, plateor frame-based monitoring of individual oysters in studies such as Garland and Kimbro (2015) and Zwerschke et al. (2016). In comparison with these methods, where oysters were attached at high densities above the seafloor on tiles or frames, the method used in this study used oysters attached to strings at a lower density of 2 oysters m −2 ; this represents a common density observed in the dredge survey During this time, individual oysters were measured for height and length using Vernier callipers. In addition, any mortalities or newly settled oysters were recorded. Strings and marker buoys were cleaned of any debris and algal growth to maintain the experiment at each monitoring event.

| Parameter estimation
Oysters monitored in the string experiment described above between March 2017 and March 2018 were used for parameter estimation (n = 177), with the initial IPM model aiming to capture the dynamics of an annual pre-reproductive census in the following year(s). In order to capture how vital rates vary with size (area, mm 2 ), a series of General Linear Models were built: a binomial error distribution (logit link) for the survival data; a Gaussian error distribution for growth increases from time t to t + 1; and Poisson error distribution (log link) for the size-based fecundity regression. Owing to the lack of data available to assess new recruit size, a mean size of 750 ± 500 (SD) mm 2 was used for new recruit size from the previous breeding season (~30 mm height). Owing to the long breeding season when new recruits may be able to grow a substantial amount between the start and the end of the season (e.g. an oyster settled in May will be larger than an oyster settled in September by the following March), a large standard deviation for recruit size was deemed necessary to incorporate the wide range of sizes of new recruits. Initial starting size distribution of the population was estimated from measurements of all oysters observed in the census dredge surveys, adjusted to the total estimated population size estimated as described above using IDW, F I G U R E 1 Schematic of the experimental design for the stringbased growth and survival experiment. The concrete blocks act as anchors, marked using buoys/fenders and oysters are directly attached to rope using Milliput adhesive. Two types of rope were used: a 5 mm polypropylene rope (rope b) was used to attach marker buoys and a 2.5 mm polyester rope used to tether the oysters (rope a) (Lown, 2019) assuming that oysters of all sizes had equal probability of being landed. To take into account the discreteness of sites within the overall population, parameters for each site were estimated individually and combined under a weighted average, weighted by the numbers present in each site in the March 2017 census. This resulted in a single IPM representing the whole BCRC MCZ area that takes account of the variation between the main sub-sites.
With O. edulis known to change sex after each breeding event, sex ratios and the percentage of the population likely to be a fecund female may vary between years. Size-based fecundity of oysters was estimated using egg counts extracted from Cole (1941) with a linear regression used to calculate height-based fecundity of ripe oysters. Owing to the potential for some oysters not to reproduce or to only reproduce as males, the percentage of fecund females within the population (i.e. females with white, grey or black 'sick') was estimated by sacrificing 149 oysters collected from the Blackwater estuary and Ray Sand within the BCRC MCZ between 22 June and 2 July 2018 to assess size-based likelihood of reproduction (mean = 73.26 ± 0.99 mm). Oysters were measured, wet-weighed and opened to check for the presence of white, grey or black 'sick' to indicate the presence of unfertilized, fertilized and developing, or near ripe eggs within the female mantle cavity (Younge, 1960). The percentage of mature females in the adult population was then included within the model as the probability of reproduction parameter or P.rep.
As with many species which include a cryptic life stage (e.g. larval stage or dormant seed stage), a common technique to overcome this 'unknown' is to implement a 'black box' called the establishment probability (P.estab). This method was used in this model owing to the lar- Here, J t+1 is the ratio of the measured population under 30 mm height at time t + 1 and T t+1 is the total census estimated population calculated through IDW. A t is the ratio of the measured popula- Owing to strings being more thoroughly checked over than the often-large dredge samples, detection of smaller sized oysters will have been reduced in census surveys despite shell substrate being available to assess.

| The IPM
The integral kernel used to create the IPM and map the size distribution at time t to the distribution at t + 1 (1 year later) was taken from (Merow et al., 2014): where z is the area in mm 2 of the oyster at time t and z′ is the area in mm 2 at time t + 1. n(z) is the size distribution of the population at time t and Ω denotes the possible range of sizes of the population. K is the full kernel, comprising P and F. P is the growth and survival kernel calculated to be: Where s(z) is the size (area of oyster) based annual survival from time t to time t + 1. g(z′| z) describes the probability density of size (z′) that an individual of size (z) can grow in a single time step, conditioned on it having survived. F(z′,z) is the fecundity kernel where: with p rep (z) the size-based probability of reproducing (assuming oysters under 1,200 mm 2 do not reproduce), f eggs (z) is the size-based fecundity, P estab the probability of an egg establishing and surviving 1 year (equation 1) and R the size distribution of 1-year-old recruits (Merow et al., 2014).
Eviction was assessed using methods described in Williams, Miller, Ellner, and Doak (2012), integrating the growth function over the bounds of the model. Initial models indicated a high probability of eviction of the largest oysters (0.3608 with error 8.6 × 10 −8 ). The maximum size was therefore increased from 1.2 times the maximum observed size (11,000 mm 2 ) to 1.8 times the maximum observed size This size range therefore remains within the natural limits of this species. Following recalibration of maximum sizes, the eviction rate of the largest oysters was <0.06 with error 6 × 10 −5 , with smaller sizes a fraction of a per cent. Therefore, no further measures were taken to account for eviction of oysters from the model.
Lambda values (λ) (i.e. the population rate of change) were calculated from the first eigenvector of K (Merow et al., 2014) and subsequent abundance estimates from model calculations were compared with total population estimates calculated from dredge survey data and IDW calculations. Confidence intervals on λ were calculated by jackknife resampling of the initial data set and recalculating λ 1,000 times.

| Model validation
To validate the model, population size distribution data were projec- To induce scenario 2 the size of new recruits was changed from 750 ± 500 mm 2 , to 2,100 ± 700 mm 2 (estimated through repeated stepwise simulations) when projecting forward from the predicted distribution of 2018 to 2019. P.estab for this year was estimated using equation 1, assuming that any oysters observed under 2,800 mm 2 were new recruits observed in the 2019 survey (with P.estab calculated as 3.63 × 10 −6 using equation 1). All other parameters remained constant.

| Parameter-based sensitivity and elasticity
In order to pinpoint which parameters are most sensitive to changes in λ, parameter sensitivity (how sensitive λ is to changes in this parameter) and elasticity (proportional sensitivity) by finite difference analysis was performed on the IPM using the initial coefficients assuming a low recruitment event (i.e. P.estab = 5.13 × 10 −7 ) and assuming high recruitment (i.e. P.estab = 8.0 × 10 −6estimated under scenario 1 above to induce high recruitment). Here, each underlying regression parameter to the IPM kernel was perturbed at random within a set threshold (0.0001) to assess how λ is influenced by changes in these parameters as presented in Griffith (2017).
Sensitivity and elasticity of λ for the model assuming low recruitment (i.e. P.estab = 5.13 × 10 −7 ) and assuming high recruitment (i.e. P. estab = 8.0 × 10 −6 ), with all other parameters the same, were calculated to determine how small changes to each parameter leads to proportional changes in λ and to investigate how these change between high and low recruitment years.

| Understanding the consequences of increased oyster recruitment: establishment probability (P.estab)
To understand impacts of recruitment, the measured abundance and with ||•|| denoting total population size calculated in each year, as also used in Metcalf et al. (2015). The population size distribution measured in the March 2017 survey, scaled to the full estimated population, was used as the starting population distribution (n0).
Confidence intervals for λ s were calculated by calculating λ s for each stochastic run and extracting the 95% confidence intervals on these values.

| Understanding the consequences of decreased oyster survival: harvesting
With fishery stock size commonly estimated using biomass as outlined above, and the current IPM model calculating population abundance, to assess fishery impacts on the European native oyster population it was necessary to convert predicted size-based population estimates to biomass estimates. Methods regarding this process and caveats associated with it are described in Supporting Information Part 2.
Using stochastic projections of the higher recruitment success Coefficients extracted from regression equations to parameterize the IPM along with establishment probability estimates calculated using equation 1 can be found in Table 1 with kernel plots found in Supporting Information Part 1.

| Parameter sensitivity
The sensitivity and elasticity of the current population is shown in Figure 4 assuming a small recruitment event in 2017. Figure 5 shows F I G U R E 5 Parameter sensitivity (above) and elasticity (proportional sensitivity) of the basic integral projection models (IPM) with zero harvesting and high establishment probability. Parameters are: surv.int (survival intercept), surv.slope (survival slope), growth.int (growth intercept), growth.sd (growth standard deviation), seed.int (fecundity intercept), seed.slope (fecundity slope), recruit.z.mean (recruit size mean), recruit.z.sd (recruit size standard deviation), establishment.prob (establishment probability, here set to 8.0 × 10 −6 ); Harv, probability of harvest of adult oystershere set to 0; p.rep, probability of being a fecund femalehere set to 0.135, i.e. 13.5% When converted to biomass, the mean log stochastic projection of the BCRC MCZ oyster populations (i.e. the stochastic model shown in Figure 6b) reaches 800 t in 30 years (range 16-66 years). Note this analysis excludes oyster populations in the four main rivers (Lown, 2019) and is based on estimations of current recruitment and survival without any management intervention.

| Understanding the consequences of decreased oyster survival: harvesting
The starting size distribution for stochastic projections incorporating harvests is shown in Figure 7  Assuming the same weight to area ratio as calculated above, a 2.5% harvest of the corresponding size distribution for adults only oysters would result in~10.8 tonnes of oysters harvested by dredge in the model starting year under current stochastic prediction with a 5% harvest indicating a catch of 21.6 tonnes, whilst ensuring λ s remains above 1. These models assume that stochastic scenario 1 is the most likely scenario, as discussed above.

| DISCUSSION
This study has incorporated extensive census data from 99 sites in a Model coefficients were found to differ between sites. Whilst no simple experiment has yet been completed using oysters from these sites, differences in growth and survival rates between areas may be due to differences in habitat, e.g. Blackwater sites have high Crepidula fornicata populations and potentially increased input of larvae from the neighbouring Tollesbury and Mersea Several Order, whereas the Ray Sand is shallow with areas of softer mud and sand and the Crouch being an area of high C. gigas density with the presence of mixed shellfish beds (Lown, 2019). There was, however, no statistical difference in temperature between sites and all sites were fully saline (Lown, 2019).

| How sensitivity of the population changes with recruitment
Parameter-based sensitivity estimates highlight the sensitivity of the current population between 2017 and 2018 to recruitment, both in high-and low-recruitment years (Figures 4 and 5;Griffith, 2017). In years when recruitment is high, whilst sensitivity of the population is still primarily driven by recruitment and adult growth, elasticity of the population (i.e. the proportional sensitivity) becomes more evenly spread across various parameters such as growth and fecundity rates. This indicates that, when establishment of recruits is high and/or regular, the population may potentially be able to withstand changes to multiple other parameters, with population rate of change influenced more evenly by a range of demographic processes as opposed to being driven by a single parameter. This has also been observed in a similar way in Californian mussels (Mytilus californianus and M. galloprovincialis), where elasticities of individual patches within a meta-population were most sensitive to juvenile and recruitment-based vital rates with individual patches most sensitive to demographic parameters than connectivity, strengthening as local retention of recruitment increased (Carson, Cook, López-duarte, & Levin, 2011;Figueira, 2009).

| Stochastic projections of recovery under interannual variation in recruitment
Stochastic projections have highlighted how the population may change under recruitment variation, assuming that vital rates such as growth and survival remain constant at the observed values of 2017.
Whilst this scenario may be unrealistic in a real-world situation, only using a single year of growth, survival data and a single year calculation for proportion of the population found to be fecund females, these projections illustrate how incorporating variable processes such as recruitment into forward projections is essential for understanding longer-term population size distributions and abundancethus managing expectations of all stakeholders. In addition, it illustrates that, when modelling species such as O. edulis, which is renowned for sporadic recruitment (Cole, 1949), conclusions require clear statements on predictive error or uncertainty for the predictions of population dynamics of future years. Historically, particularly large recruitment events in this species have occurred once every 7 years or less (Cole, 1949). Stochastic run 1 assumes no variation in growth and survival and incorporates estimated recruitment rates that assume a medium-sized recruitment event in 2017. This recruitment F I G U R E 8 Forward stochastic projections projecting forward from the extracted 'recovered' population selecting establishment probability at random from 5.13 × 10 −7 , 8.0 × 10 −6 , 3.79 × 10 −7 and 2.22 × 10 −6 and inducing 0, 1, 2.5, 5 and 10% probability of mortality of landable sized oysters every year. The red line highlights mean log population size for all model runs with black lines representing individual model runs subsequently re-occurs semi-sporadically at random, assuming that Whilst it is recognized that higher levels of recruitment have been induced to replicate an observed distribution of potentially 2-year-old recruits (i.e. scenario 1), higher levels of recruitment than those induced here have previously been observed in the Essex estuaries. A previous stock survey of the River Crouch showed 80% of the population as 1-year-old recruits (Shelbourne, 1957 Conversely, other previous studies have indicated no recruitment in multiple years, highlighting that the reverse may also be true (Bromley, McGonigle, Ashton, & Roberts, 2015). Owing to climate warming, oyster growers in Essex now expect reasonable recruitment events to occur at least once a year, with particularly large recruitment events occurring on intermittent occasions (Essex oyster growers, KEIFCA stakeholder workshop July 2019), as has also been observed in scallops (Shephard, Beukers-Stewart, Hiddink, Brand, & Kaiser, 2010). In Essex, the grounds of the private oyster-growing areas are specifically managed to promote recruitment and it is therefore reasonable to expect recruitment rates within the public grounds (i.e. the grounds outside of the several orders and private fisheries) to be lower than those observed in the several order and private fishery grounds. Likewise, from this comparison with private managed areas, it is also reasonable to assume that management to improve recruitment success is feasible in the wider MCZ. This may include desilting or spatting substrate habitat improvement activities (Laing, Walker, & Areal, 2005).
Understanding cryptic life stages is highly problematic for population modelling. Cryptic life stages are found in species ranging from fin fish to plants and include any species with a larval stage, small seed dispersal phase or seed bank. This is commonly overcome with the use of a 'black box' , described here as the probability of recruitment (Merow et al., 2014 Ostrea edulis are notorious for highly fluctuating populations and a succession of chance low or high recruitment events will greatly alter the population size and extent, with fisheries often able to be sustained by one or two particularly large recruitment events for multiple years (Spärck, 1949). Historical oyster populations have taken in excess of 20 years to be deemed adequately restored for fishing to occur with several years of highly favourable recruitment required for populations to recover (Spärck, 1949). The estimate of 30 years to recover from~300 to~800 t may therefore seem appropriate given the extent of population growth required to open a sustainable fishery, noting that the increased frequency of warm summer months means that southern oyster populations are now experiencing spawning conditions in most years. It is possible that contemporary warming (up to some limit), a super-recruitment event and any active intervention could result in faster recovery times than are estimated here.
Whilst making predictions on the number of individuals is key for assessing the recovery of a population, understanding how other ecosystem services may benefit is also necessary to determine the long-term viability and benefit of a population. As with other species of oyster, O. edulis shells show higher species diversity than other non-living hard substrate (Smyth & Roberts, 2010). In addition, an increasing density of native oysters results in greater species richness, in the absence of the invasive slipper limpet Crepidula fornicata (Lown, 2019). These metrics, such as biodiversity benefits, are essential in fully determining a recovered ecosystem as a whole.

| Modelling harvest scenarios
Once recovered (i.e. once the population reaches 800 t or~10 million oysters), λ s values obtained inducing various levels of harvesting of adult populations indicate that a population of this size would result in~3 million oysters of landable size. A 2.5% harvest or less of this population would result in~75,000 oysters of a range of sizes able to be landed and would enable harvesting to occur whilst maintaining λ s > 1. A 5% harvest would maintain λ s close to 1 (1.0018) and result in 150,000 oysters. This equates to~10 t of oysters per year for a 2.5% harvest and 20 t for a 5% harvest.
Assuming that the maximum level of 250 kg native oysters is harvested on any one day, this would allow 40 boats to be granted a single day licence for the fishery for a 2.5% harvest rate (or 10 boats for 4 days assuming the full catch allowance is landed per day) (KEIFCA, 2019a). This is under the assumption that growth and survival rates do not decrease (or increase) and establishment of new recruits is not below an average of 2.53 × 10 −6 over 4 years.
This highlights that, if it is maintained at low levels within a growing, healthy population with regular recruitment, fishing within sustainable guidelines may be possible, provided there is sufficient habitat available to sustain the population. Further research, particularly in understanding habitat limitations, density dependence and drivers influencing the establishment probability of new recruits, is required in addition to increased years of data collection to further validate our model and improve its predictive power.

| CONCLUSION
If current estimates of growth, survival, fecundity and recruitment are accurate and do not improve further in the absence of active intervention recovery of the BCRC MCZ, the native oyster population is estimated to recover in 16-66 years with a mean of 30 years. Following recovery, this study has shown that it is reasonable for a sustainable fishery to be developed with harvest rates <5% that would meet local stakeholder needs such as sufficient catch allowance to justify a small fishery, whilst ensuring sus- essexnativeoyster.com). Owing to staffing and resource changes, 2019 data collection was interrupted. Owing to the coronavirus pandemic, 2020 data collection has been interruptedthanks are due to all our partners scrambling to maintain data collection in difficult conditions.