European native oysters and associated species richness in the presence of non‐native species in a southern North Sea estuary complex

There are growing calls to restore populations of European native oysters (Ostrea edulis), on the premise that restored populations will support a range of ecosystem services with an emphasis placed on restored oyster habitats promoting biological diversity, however benefits associated with naturally occurring O. edulis remain unclear. We undertook biannual surveys in the Blackwater, Crouch, Roach and Colne Estuaries Marine Conservation Zone (BCRC.MCZ), a highly sedimented estuary complex in the southern North Sea, to investigate links between natural densities of O. edulis (0–4.2 m−2), the prevalence of other dominant habitat features such as non‐native slipper limpet (Crepidula fornicata), dead shell abundance and epibenthic macroinvertebrate species richness. Increased epibenthic species richness was associated with O. edulis, even at densities below the OSPAR Commission recognized definition of an oyster bed (5 oysters m−2). Our analysis predicts increased associated species richness with density of native oysters (e.g., +1.6 additional species at 1 oyster m−2 or + 2.8 species at 5 oysters m−2), but only in areas with lower density of C. fornicata. Where C. fornicata are at higher density, the potential benefits of oyster restoration for associated species were curtailed. This may explain the observed asymptotic relationship between oyster density and diversity at 1 oyster m−2. In these and other high Crepidula density areas we recommend extending native oyster habitat even at low density. This may be of particular interest to areas with the protozoan oyster parasite Bonamia ostreae, which spreads more easily in high‐density areas. These lower density thresholds should also be considered for future management decisions—closing harvests so they do not reduce density further and impair biodiversity services of the habitats. In conclusion, while C. fornicata may be a useful oyster settlement substrate, we find that it limits the potential increases in associated species gains of oyster restoration.


| INTRODUCTION
The past two decades have seen increasing attention toward rewilding and restoring habitats with efforts to re-establish extensive natural processes that drive ecosystem functions (Pereia & Navarro, 2015). Contemporary restoration is largely based on the creation or enhancement of protected areas, often directed toward the reintroduction of ecosystem engineers, with an emphasis on terrestrial and freshwater ecosystems (Law, Gaywood, Jones, Ramsay, & Willby, 2017;Ockendon et al., 2018). Marine-based restoration projects have largely focused around tropical and sub-tropical corals and seagrasses, and restoration of temperate marine benthic habitats is now gaining momentum on grasses, marshes, and mudflats (Airoldi & Beck, 2007). Despite this, the need to restore more structurally complex habitats has been noted (Airoldi & Beck, 2007). Coastal restoration projects are increasingly common, with oyster restoration an emerging field, specifically in the USA, Canada and Australia, where environmental, economic and cultural significances are being recognized (Coen et al., 2007;McLeod, zu Ermgassen, Gilles, Hancock, & Humphries, 2019).
Globally, oysters and their habitats are estimated to have declined by 85% (Beck et al., 2011) with long-term declines in the European native oyster, Ostrea edulis, reported since the industrial revolution (Spencer, 1990). O. edulis has supported fisheries across Europe since Roman times, remaining a culturally important marine species to this day, however, has only recently become a focus for extensive marine habitat restoration. The species is listed as a Priority Species in the Post-2010 UK biodiversity framework (Pogoda, Brown, Hancock, & von Nordheim, 2017;Syvret & Woolmer, 2015). These include calls for widespread restoration and associated socioeconomic benefits like fishing and cultural heritage (Fariñas-Franco et al., 2018;Helmer et al., 2019;Pogoda, 2019).
The focus of our study addresses another common justification for oyster restoration-the recovery of associated biodiversity. The predominant literature on benefits of oyster restoration to biodiversity are based on faster growing, multi-dimensional-reef forming species such as the eastern oyster (Bersoza Hernández et al., 2018), with the benefits often unquantified for the majority of other species of interest (zu Ermgassen et al., 2020). Shellfish as a group are known to enhance marine biodiversity by creating hard substrata (Smyth & Roberts, 2010), however, species-and scale-specific benefits of restoration need to be addressed (Gillies, Crawford, & Hancock, 2017;Pogoda, 2019;zu Ermgassen et al., 2020). It has been demonstrated that live O. edulis shells can host more epibiont richness than other hard and nonlive surfaces (Smyth & Roberts, 2010), however, this will not necessarily scale up to increased marine biodiversity at larger spatial scales or with a broader range of taxa. Experimental approaches have assessed species associations of Ostrea and Crassostrea species, as well as differences between them in enclosed structures, cages and mixed species beds with no observed differences in assemblages between manipulated O. edulis versus C. gigas densities (Zwerschke, Emmerson, Roberts, & O'Connor, 2016) [Correction added on 16 January 2021, after first online publication: the word "species" before "Ostrea" was removed in the sentence "Experimental approaches have assessed species associations of species Ostrea and Crassostrea species.."]. This points to the shellfish substrate or reef structures bringing real benefit to biodiversity rather than species specificity per se, with an estimated 60% increase in species richness observed in mixed shellfish beds over smaller spatial scales (Christianen et al., 2018). Recent studies of historical catch records and literature suggest such larger scale species richness associations with native oysters did once exist (Bennema, Engelhard, & Lindeboom, 2020). Therefore, it should be expected, but there remains to be an assessment of whether diversity associations form with naturally occurring subtidal Ostrea oysters at larger spatial scales (Gillies et al., 2017;Mcleod et al., 2019). The knowledge gap is particularly pertinent across the soft subtidal marine and estuarine sediments found across Europe where European native oysters were once abundant (Bennema et al., 2020), as our understanding of large scale oyster-biodiversity associations from a range of other species are from intertidal habitats (Green & Crowe, 2014). Understanding where density thresholds for biodiversity associations may occur is also important in the context of forming management objectives about restoration (Guy, Smyth, & Roberts, 2019). If project aims are to boost biodiversity while supporting socioeconomic recovery of Ostrea oyster fisheries, or to increase sustainability of existing fisheries, adoption of density thresholds may be appropriate to restrict fishing when oyster densities approach those below which biodiversity may decline (KEIFCA, 2018;Lown, Hepburn, Dyer, & Cameron, 2020;McLeod et al., 2019).
An important confounding factor is the presence of the non-native slipper limpet (Crepidula fornicata) in many European estuaries where oysters are present. Slipper limpets inhabit areas once classed as subtidal mixed sediments and native oyster beds, reducing the likelihood of natural recovery of Ostrea populations (as they would have been 100-1,000 years ago) (Blanchard, 1997;Helmer et al., 2019). It is therefore important to understand how the presence and abundance of C. fornicata affects the relationship between native oysters and associated species.
To address these knowledge gaps, we present the first study on diversity associations across naturally occurring densities of O. edulis over an extensive protected area, spanning 28,400 ha, designated to protect and restore O. edulis, their habitat and associated heritage fisheries. Using data from multi-year and multi-season sampling over a range of densities, including some of the highest densities on the UK coastline, we address how oyster density, and other sources of hard substrate such as slipper limpet or dead shell are associated with macrofaunal species richness variation.
We test the hypothesis that O. edulis are more strongly associated with increased species richness than that provided by other hard shell or non-native shellfish species alone. This study will help to inform national and international efforts toward Ostrea oyster restoration, helping to define what good ecological status of oyster habitats looks like and develop policy and best-practice, specifically regarding interpretation of oyster density thresholds for management decisions (Haelters & Kerckhof, 2009).

| MATERIALS AND METHODS
The Blackwater, Crouch, Roach, and Colne Estuaries Marine Conservation Zone (BCRC.MCZ), Essex UK was selected for this study following the designation of this protected area for the presence of native oysters and native oyster beds (for further site information see Supplementary information Part 2).
Initial surveys were performed in September 2014 and 2015 to assess the distribution of O. edulis and C. gigas within the MCZ. In initial surveys, the whole BCRC.MCZ was split into grids measuring 0.5 0 latitude × 1 0 longitude (0.925 × 1.85 km) with each grid square sampled at the center point (Wiggins, 2014). If native oysters were present, four additional "sub-site" dredges were completed in the Northwest, Northeast, Southeast and Southwest of the grid, resulting in five dredge tows of 120 m 2 each per grid. Subsequently, we repeated surveys during winter (March/April) and summer (September/ October) in 2016 and 2017 to assess oyster abundance, distributions, and relationships with subtidal epibenthic associated species. These repeated surveys in 2016 and2017 only included grid cells and associated sub-sites where oysters of any species were found in 2014/15 (n = 103 sites; Table S2). Only these latter surveys are used in our analysis. Seven sites were not surveyed in winter 2016 due to boat repairs but were sampled in subsequent surveys. All sites were fully saline and shallow, with site depths ranging between 0 and 16 m below chart datum. Between June 2017 and June 2018 surface temperatures were recorded at the three main oyster bed areas (Blackwater, Crouch and Ray sand). Subsequently, benthic Hobo loggers were used to record daily midday temperatures between June and September 2018. Previous work, and our own has found no significant effect of water temperatures on oyster aggregations at this scale (Allison, Hardy, Hayward, Cameron, & Underwood, 2020).
All survey samples were completed with 1.2 m wide ladder dredges (Table S3) towed for 100 m at a ground speed of 2 knots in the center point of each grid square or sub-site location (Wiggins, 2014). Dredge contents were sorted into constituent parts including dead shell and live C. fornicata, O. edulis, and epibenthic species. We identified and recorded abundance of epibenthic species present to as low a taxonomic resolution as possible in the field (Gerwing et al., 2020), including abundance and weight of native oysters and weight of live C. fornicata (Table S2) (See Table S4 for full species list).
Dredge use was necessary to assess the abundance of O. edulis over an extensive area and provide measures of oyster density to develop restoration and conservation objectives. Alternative sampling was explored, for example, grab sampling, but found to not be an accurate representation of densities of subtidal epibenthic macrofauna, including shellfish, compared with dredging (e.g., Cameron, 2017). All species (except C. gigas) were returned to the sea alive close to their collection site.

| Scaling of the predictors of species richness
The objective of statistical models of species richness were to predict the change in taxa diversity with variation in abundance of different hard substrata. Such a statistical modeling approach is required due to the highly correlated nature of the habitat explanatory variables, where a univariate or basic regression approach would be inappropriate. In all models, native oyster counts and two measures of the availability of other hard substrata, total dead shell weight, and live C. fornicata weight were used as explanatory variables (i.e., predictors), with species count data as the total number of distinct taxa from the species list within the dredge. Fish, alga, and some other epibionts were not included as were either not suitably quantified by dredge sampling or could not be processed at sea. As it was not possible to quantify the patchy and variable densities of native oysters over small distances, predictors were converted to average densities per m 2 across the 120 m 2 dredge.
It is recognized that dredges and trawls only sample a fraction of the fauna on the seabed (Eleftheriou and Moore 2013), therefore a review of dredge efficiency literature between Natural England, Kent and Essex Inshore Fisheries Conservation Authority (KEIFCA) and ourselves led to efficiency estimates. Average native oyster abundance was adjusted to a 20% dredge efficiency, where it is assumed 20% of the native oysters are caught in the first pass (Allison, 2017;Dijkema, 1983;KEFICA, 2016). This 20% efficiency was also applied to live C. fornicata (kg/m 2 ) with stacks of C. fornicata being a similar weight and size to O. edulis. Dredge efficiency has previously been calculated to be approximately twice as efficient for capturing live eastern oysters as it is for cultch (Taylor et al., 2007). Whilst it is recognized these two species of oyster occupy different niches and grow at different rates, eastern oysters remain a similar size and weight as O. edulis, therefore a dredge efficiency of 10% was used for dead shell average density (kg/m 2 ). All future abundances for predictors are stated as scaled up average densities incorporating these efficiencies. Calculating dredge efficiency is notoriously difficult, therefore a corresponding analysis of un-scaled predictor variables can be found in Supplementary Information 7. Scaling up did not impact overall predicted trends and statistical significance, however, it did alter the value of the coefficients ( Figure 1).
Collinearity between predictors was assessed using Variance Inflation Factors (VIF) (Zuur et al. 2013). All VIF statistics were below 1.5 therefore all predictors were retained for regression analyses ( Figure S2). Rarely sampled taxonomic groups such as annelids and small decapods (e.g., Crangon crangon) were grouped in terms of presence/absence and termed "other species" for individual species analyses. Species such as Sabellaria spp. and barnacles were only measured as presence/absence due to inaccuracies in counting individuals, with Sabellaria largely broken up in the dredge and barnacles often found in extremely high numbers (for full taxa list and usage see Table S4).

| Model selection
All statistical analysis was completed in R version 3.4.3 (R studio version 1.0.143) (R Core Team, 2017). To visualize the effects of oyster density on epibenthic biodiversity, number of taxa observed were plotted as a function of average O. edulis density (m −2 ) per dredge with a loess smoother and 95% confidence interval. This was repeated for the other explanatory variables-total shell (i.e., all shell caught in the dredge) (kg/m 2 ) and total weight of F I G U R E 1 Map of all sites surveyed in the UK native oyster survey with Kent and Essex Inshore Fisheries Conservation Authority. All start points of dredges from 2014 to 2017 have been plotted along with the Blackwater Roach Crouch and Colne Marine Conservation Zone and Several Order/Private fishery boundaries. Many dredges occurred at the same point resulting in an overlap of map coordinates. The most recent survey where the area was sampled is on the top layer of the map live C. fornicata (kg/m 2 ). Due to the highly correlated occurrence data for both explanatory and response variables-these trends or univariate approaches must be interpreted with caution and a modeling approach that can account for the co-occurrence is required.
Variation in (i) number of associated taxa for all seasons, (ii) number of associated taxa for summer, (iii) number of associated taxa for winter and (iv) taxa presence/absence composition, as a function of the three substrate explanatory variables were estimated using multivariate generalized linear models using statistical package "mvabund" (Wang, Naumann, Wright, & Warton, 2012). This gave a community level metric not constrained by the distribution of the data (ManyGLM, i, ii, and iii: negative binomial response and iv: binomial response for the taxa richness and presence, respectively, Wang et al., 2012). Model distribution fit was assessed with plotted Dunn-Smyth residuals. All models used 999 bootstrap permutations. In winter 2016, some sites could not be sampled resulting in unequal number of samples taken in each survey. Subsequently "case" resampling was used to incorporate repeated measures designs of multiple surveys with differing number of samples (Wang et al., 2012). All predictors were added to the model incorporating a fully factorial design to identify potential interaction effects between different hard substrata types. Initial analysis included C. gigas however, as the number of dredges containing C. gigas were rare compared with other predictors (n = 126/397) and, following nonsignificance with species diversity in preliminary analysis, this species was removed as a predictor and included in the associated species number defined above. Starting in all cases with an apriori maximal model of all main effects and their interactions, the best minimal model was selected with stepwise AIC where all two-way interactions between predictors were retained within Models i, ii, and iii, and main effects only retained for Model iv. Residual plots can be found in Supplementary Information 9. Coefficients were extracted from best models to determine the impact of predictors on individual species (see Supplementary Information Part 17 for full R code). Respective predictions were then made on a dummy data set using Model i (defined above), assessing variation in the number of associated taxa for all seasons.
Uncentered counts of native oysters were used in all analysis. To test the robustness of our results, a comparative analysis of each of the models discussed above was repeated using all predictors standardized to kg/m 2 and for centered predictors (e.g., standardized predictor variables to clarify interpretation of the intercept). Comparative analysis can be found in Supplementary Information 11. These variations did not alter our conclusions. (c) Weight of live Crepidula fornicata (kg/m 2 ) (adjusted to 20% dredge efficiency) versus number of different species within the dredge. All plots show a loess smoother with span = 1 and 95% CI have been added for ease in observation. Each of oyster, dead shell and C. fornicata density are highly correlated and require a different approach to elucidate which are the strongest predictors of associated species richness

| RESULTS
A total of 396 dredges were completed between February 2016 and October 2017, identifying 39 taxa to be included in our analysis (Table S4). Oyster sizes ranged from 6 to 119 mm shell height (average of 66.25 mm ± 0.295 SE, for size distribution of oysters see Figure S2).
3.1 | Does oyster, Crepidula, or shell density influence number of observed species?
Increases in the abundance of O. edulis density, C. fornicata density and dead shell weight were associated with significant increases in the number of taxa present, from an average of three taxa in the absence of each environmental variable, to an asymptote at 5-7 species (Figure 2). Most deviance was explained by dead shell weight (LR = 75.87, p < .001), then C. fornicata (LR = 20.44, p < .001) and the native oyster density (LR = 12.03, p < .001). Significant interactions between average dead shell and C. fornicata weight (LR = 10.91, p = .003), dead shell and O. edulis average density (LR = 15.24, p < .001) and C. fornicata weight and O. edulis density (LR = 5.14, p < .001) were observed. Coefficients extracted from these models were plotted to identify the impact of an increase of 1 kg of any substratedead shell, live C. fornicata, or one additional O. edulis m −2 on species richness. Increases in observed associated taxon number were driven by increased volumes or density of substrate types in isolation while interactions between substrate types and their abundance resulted in negative impacts on the number of associated species (Figure 3; Supplementary Information 13 and 14 for analysis using weights of oysters as predictors as opposed to counts and a comparative analysis using mixed models. Alternative models made no qualitative difference to our conclusions).
To visualize the relative importance of the effects of biodiversity predictors, a predictor dataset using incremental increases in O. edulis density, C. fornicata weight, and dead shell weight were used to predict changes in number of associated taxa m −2 using the original model ( Figure 4); for example, the number of associated species observed with increasing dead shell or C. fornicata abundance in the absence of native oyster. While total shell explained most deviance, it has a relatively small effect on biodiversity compared with other variables (Figures 3  and 4). This is reflected in the model predictions (Figure 4). A key prediction is that consistently large increases in associated species diversity were observed with increasing average O. edulis abundance in the absence of C. fornicata (e.g., +1.6 additional species at 1 oyster m −2 or + 2.8 additional species at 5 oysters m −2 ). However, in the presence of C. fornicata, the positive effects of oysters on biodiversity were reduced (Figure 4).
The same predictors were significant in both summer and winter, with coefficients largely the same as the annual analysis and did not change between seasons (for more information see Supplementary Information 14). However, changes in the number of taxa observed indicate increased diversity in summer surveys (LR = 5.303, p = .021). A simplified analysis that discretizes oyster density into bins and treats it as a factor, to help clarify where the greatest changes in biodiversity occur are found in Supplementary Information 15, however, this did not qualitatively affect the results.

| Does oyster density influence species composition?
In the analysis of species composition (i.e., Model iv), average dead shell, live C. fornicata weight and O. edulis density were all found to have a significant effects on the composition of all species. (Dead shell (kg/m 2 ) F I G U R E 3 Coefficients with SE extracted from negative binomial ManyGLM identifying effect of total shell weight (kg/m 2 ), live Crepidula fornicata weight (kg/m 2 ), O. edulis density (m −2 ) and interactions between total shell: live C. fornicata weight, total shell: O. edulis density and also C. fornicata weight: O. edulis density with intercept value on the total number of species observed in a dredge. Respective dredge efficiency percentages used are shown in brackets with densities calculated over an average of 120 m 2 . NB different units for effect size of O. edulis on species richness-currently 1 kg of C. fornicata gives the same biodiversity benefit as a single O. edulis m 2 LR = 196.0 p = .045, Live C. fornicata (kg/m 2 ) LR = 106.5, p = .039 and O. edulis density (m −2 ) LR = 105.7, p = .021). There were no significant interactions therefore these were removed. Unstandardized coefficients for the model for individual species occurrence response are plotted in Figure S9, highlighting increases in density of anemone species, C. gigas, queen scallops, brown and velvet crabs and clam species are more likely to be present with increasing native oyster abundance.

| DISCUSSION
O. edulis may be a species that could boost biodiversity through acting as an ecosystem engineer, should it recover its lost range and abundance. This is due to hard shellfish being a preferable substrate to epifaunal species (Smyth & Roberts, 2010). The majority of studies to date have, however, assessed the associations of other oyster species (Bersoza Hernández et al., 2018), occurred before the prevalence of certain diseases and invasive species (Mistakidis, 1951), or used manipulated densities in unnatural/manipulated settings (Zwerschke et al., 2016).
Here we have shown that whilst the presence of O. edulis is associated with increased diversity up to an observed average density of 1 oyster m −2 (24 oysters per experimental dredge, Figure 2), no further increases are predicted across the majority of the habitats we sampled (with average densities up to 4.2 oysters m −2 and 1.5 kg slipper limpets in a 120 m 2 dredge). Greatest increases in species diversity in situ were found between oyster-free areas and areas of low average density (0.5-1 oyster m −2 ), indicating significant increases in associated species are found at oyster densities below those required to classify an oyster bed by the OSPAR definition of 5 oysters m −2 (Haelters & Kerckhof, 2009).
Whilst our analysis demonstrates an additional 1.5 species at 1 oyster m −2 and predicts an additional 2.8 species at 5 oysters m −2 , for comparison with other literature, this still represents an 87% rise in associated species-similar to biodiversity associations with mixed native and non-native shellfish beds which noted a 60% rise in diversity (Christianen et al., 2018). It should be noted that we have only sampled a subset of potential species diversity in the current study (e.g., macrofauna invertebrates). In particular, we found anemone species, F I G U R E 4 Predictions in changes in number of species per m −2 with standard errors extracted from the negative binomial ManyGLM model. Predictors are number of native oyster m −2 , weight of Crepidula fornicata kg/m 2 and weight of total shell kg/m 2 . Columns of plots are split by incremental weight of C. fornicata (0-1.5 kg/m 2 ) with rows split by incremental total shell weight (0-2 kg/m 2 ) C. gigas, queen scallops, brown and velvet crabs and clam species are most likely to be associated with increased native oyster abundance ( Figure S9) suggesting potential socio-ecological benefits could arise from recovery of extensive and complex native oyster subtidal habitats. Assessing the presence and prevalence of sponges, bryophytes, algae, and so forth will improve on future predictions. In addition to including more taxonomic groups, for example, algae or fishes, it is possible a low regional species pool or dredge selectivity affects the maximum associated species richness in our study and so a range of sampling methods could be employed in any future assessments.

| Potential for density dependent associations with biodiversity
Notably, further increases in associated species above 1 oyster m −2 /100 m dredge in the raw MCZ data were rare. This could be due to a variety of reasons. Firstly, increases in associated species with higher native oyster densities in raw plots are not observed due to negative effects caused by high densities of C. fornicata. Secondly as areas of higher oyster density are relatively rare, while they are common in comparison to other sites in the UK (Clark, 2017), this may be limiting patterns of higher associated species richness emerging in our study. Finally, O. edulis may be found in habitats already supporting increased biodiversity due to high food availability, hydrodynamics, or some other unmeasured variable, without directly influencing diversity of species.
Our predicted estimates of the number of additional species under varying scenarios show a linear increase with increasing native oyster density at low C. fornicata density, however, we found these increases are suppressed and even reversed when C. fornicata are at high density (Figure 4). Colonisation by C. fornicata results in a homogeneous seafloor which alters the community type (Blanchard, 2009). This has also been observed in Mytilus edulis and C. gigas beds, with decreased species richness observed at highest density samples (Beadman, Kaiser, Galanidi, Shucksmith, & Willows, 2004;Green & Crowe, 2014). Reductions were attributed to competition for space, physically covering the sea floor resulting in a reduction of predator species which feed on infauna due to inaccessibility to soft sediment, and potential variation in hydrodynamics caused by the change in shell substrates (Green & Crowe, 2014). This may be occurring in the Essex estuaries however, percentage cover was not evaluated in this study. We would note that in a soft mud dominant estuary system there may still be biodiversity benefits, however small, of the hard surface mosaic that slipper limpets can create (De Montaudouin & Sauriau, 1999). We have shown this positive effect of slipper limpets on biodiversity in our model predictions in the absence of O. edulis however this effect is reduced in comparison to positive effects observed with O. edulis or shell, for example, from 0 to 1 oyster (e.g., max 100 g) results in +0.3 species in the absence of other hard substrate, however, from 0 to 0.75 kg slipper limpet results in a prediction of +0.25 species ( Figure 4). As with many benthic species, slipper limpets can re-emerge following burial by light sedimentation (Powell-Jennings & Callaway, 2018). The BCRC.MCZ is a muddy estuary with concentrations of over 50 mg/l 1 suspended sediment (Moffat, 1995). Nonalgal suspended solids above 10 mg/l 1 are considered normal for the estuaries of other major European rivers but the Thames and Anglian plume are known for values above 30 mg/l 1 (CEFAS, 2016). Given their abundance it is therefore likely C. fornicata are able to provide a constant influx of hard substrate, providing clean shell for settlement of new species year-round, including European native oysters. The context of where C. fornicata is a pioneer species for biodiversity, including in management for the recovery of native oyster, or where they prevent maximum biodiversity gains from oyster restoration requires further research. Such management is possible as Essex native oyster producers regularly state that slipper limpets can be useful to establish suitable oyster settlement habitat, if that habitat is managed and prevented from becoming homogeneous.
Due to the rarity of high density O. edulis areas (8 dredge samples from a total of 396 have an average density over 2.5 oysters m −2 over the 100 m dredge tow), to validate conclusions taken from our analyses, we repeated our dredge survey with the same methods in high oyster density areas in a managed private oyster several order near our BCRC.MCZ study area (Allison et al. 2020). These areas had similar oyster density ranges and associated species numbers (3-12 species with a range of 0.54-4.6 oysters m −2 within the private aquaculture areas and 0-15 species with a range of 0-4.1 oysters m −2 within the wider BCRC.MCZ). This suggests that high biodiversity-high oyster density relationships may exist but something, probably slipper limpets, limits their occurrence. This hypothesis would benefit from being tested experimentally.
While ideas presented above for limited biodiversity gains in mixed native oyster and C. fornicata beds remain hypotheses that need to be tested, our analysis predicts biodiversity gains of restored oysters when C. fornicata abundances are low. Therefore, reductions in C. fornicata may be an appropriate management tool to explore to maximize biodiversity associations from European native oyster restoration. We recommend that a restoration and fishery management objective should be to restore and maintain European native oysters of at least an average density of 1 m −2 across a 120 m 2 dredge. This is below the oyster bed definition of 5 oysters m −2 , however, should improve associated species diversity while taking into account the cautions raised above such as disease risk with diseases such as Bonamia spread more easily in high density oyster beds (Doonan, Cranfield, & Michael, 1999). To maximize marine biodiversity recovery, if sitespecific funding or broodstock availability were to be an issue, it may be more strategic for restoration projects to maximize the footprint of a project at low to moderate oyster density, boosting the suitability of extensive habitat for oysters and associated species instead of smaller areas with higher densities.

| Marine sampling challenges and observing real density dependent biodiversity relationships
While we observed average native oyster densities up to 4.2 oysters m −2 , like all dredge sampling we cannot measure absolute density in every square meter. However, we assume native oysters are likely to have a patchy nonuniform distribution and, subsequently, there will be areas, which surpass the OSPAR definition of 5 oysters m −2 in the BCRC.MCZ. Despite these challenges, we would strongly emphasize a focus on both optimism and caution when presenting the case of O. edulis restoration. We suggest optimism in promoting moderate levels of species richness as an achievable goal for restoration, as increases in species richness can occur at relatively low oyster densities within safe disease-risk limits (up to 1.26 oyster m −2 , Doonan et al., 1999). Some of those associated species such as we found in this study, Brown crab, support the economy of inshore coastal communities. However, we also advise caution in promising that restored natural densities of European native oyster will deliver substantial increases in subtidal benthic diversity. Protection of high-density oysters as "beds" could create unnecessary conflict due to policy implications of native oysters as a species versus a habitat regarding fishery management. This conflict can be minimized through adaptive and spatial management, and consideration that ecological restoration and fishery "recovery" are two different but complementary achievable objectives (Lown et al., 2020). Our study does not yet incorporate other ecosystem services (e.g., denitrification or fish nursery or foraging habitat potential) or minimum densities for successful reproduction, which may also affect the oyster density thresholds at which benefits to society accrue.

| CONCLUSIONS
This study has demonstrated a positive association between the presence and density of O. edulis and an increased diversity of sublittoral benthic communities. This has implications for the management of this species, in terms of restoration and any future harvesting, particularly in areas where C. fornicata are prevalent. We cannot infer why this positive association with O. edulis occurs. Additional biodiversity benefits associated with highest density oyster areas were rare in the current conditions of the BCRC.MCZ. We have shown this is likely to be driven by high C. fornicata abundances in many sites, predicting steeper biodiversity gains from oyster restoration in low C. fornicata density areas. We recommend exploring outcomes of oyster restoration in areas of varying C. fornicata density to inform management and practice, and determine whether density reductions of C. fornicata prior to restoration efforts are worthwhile. Under current conditions of high C. fornicata abundance and acknowledging this is a recommendation that ignores other potential motivations for restoration, increasing the area over which oysters inhabit is likely to be more beneficial to epimacrofaunal sublittoral biodiversity than boosting densities over smaller areas in order to meet what may be an arbitrary oyster bed definition of 5 oysters m −2 .
DATA AVAILABILITY STATEMENT Data are available in the University of Essex Research Repository associated with this publication title. No exact location data will be associated with any of the species data provided. This is due to the protections associated with current stocks of Ostrea edulis within our study area.

ETHICS STATEMENT
All work was undertaken with dispensation to work within the BCRC MCZ from all appropriate and competent authorities. None of the research on animals in this publication required UK Home Office approval as was below ASPA thresholds.