Picking a way forward: valuing and managing traditional shell ﬁ sh gathering for Littorina littorea

– Littorina littorea (periwinkles) have been harvested by hand picking from the shore since prehistoric times. Harvests are generally unregulated, catches are not linked to particular shores and ﬁ sheries statistics are considered to be unreliable. The absence of key data has made it dif ﬁ cult to develop harvesting recommendations. Surveys around Strangford Lough, Northern Ireland were used to investigate the size structure and relationships among densities in different size classes. Three size classes were identi ﬁ ed in surveyed L. littorea , with mean shell lengths of 0.81, 1.56 and 2.48 cm. Assuming that the age classes represent year classes, data across different shores suggested that the ratio between densities in successive year classes was not constant. Proportionally fewer individuals were found in the larger, older, size class as the density of the smaller size class on a shore increased. This density-dependent relationship was modelled with a Ricker curve for the year 1 to year 2 and the year 2 to year 3 transitions. The predicted transition rates from Ricker curves were used in a size-structured model to describe L. littorea dynamics. An emergent property of the size-structured model is a decline in mean shell length with overall density of a population. This prediction was supported by the survey data from Strangford Lough and by an independent survey of Irish shores. The size-structured model predicts potential harvests of individuals above 2.06 cm as a function of recruitment rate. Maximum harvest was predicted for a density of 5 year 1 individuals m (cid:1) 2 , leading to 13.8 year 3 individuals m (cid:1) 2 or an estimated annual harvest weight of 67 g m (cid:1) 2 . Modelled estimates of production provide a means to value shores and develop harvest predictions for management purposes.


Introduction
Littorina littorea (L.), the common periwinkle, is a conspicuous intertidal snail found in the North Atlantic. Periwinkles have a long history of exploitation, shells are commonly found in prehistoric shell middens (Gutieŕrez-Zugasti, 2011;Hood and Melsaether 2016). When harvested, L. littorea are still mostly gathered by hand from the shore, usually for food, but occasionally as fishing bait (Fowler, 1999). Landings figures are generally considered unreliable (Tully, 2017) and to underestimate the true catch (Chapman, 2006). Recent landings figures from Ireland are typical in showing large fluctuations between years: mean 808 (SD 573) tonnes per year between 2004 and 2015 (Tully, 2017). Based on a sales note (first sale) price of €2.04 kg À1 (Anon, 2014) the L. littorea catch in Ireland has an annual average value of approximately €1.6 million. Value of the UK catch is likely to be of a similar magnitude. Harvesting is an informal activity used to supplement income and may therefore be affected by the availability of other sources of earnings and wider economic conditions (Cummins et al., 2002;Chapman 2006).
A number of factors combine to make the harvesting of periwinkles a data-poor fishery, stemming from the informal nature of the harvest. The boundaries of stocks being exploited are difficult to identify, as detailed records are generally not available for harvest locations. Locationspecific trends in population abundance and effort are generally unknown. Harvesters pick selectively and typically fill bags that are collected from a number of shores before the catch is taken to wholesalers. The lifecycle of L. littorea has a planktonic larval stage. Reproduction involves internal fertilization with egg capsules released into the plankton. Following hatching, the larvae spend several weeks in the plankton and the average dispersal distance has been estimated to be 22.9 km from a reanalysis of allozyme data (Kinlan and Gaines, 2003). The link between stock and recruits is therefore complex and likely to vary with location as the supply of new individuals to an exploited population could originate from different source populations.
As grazers, L. littorea can influence the ecology of shores. Lubchenco (1983) described preferential consumption of ephemeral green and red algae by L. littorea that facilitates the establishment of fucoid algae. The recruitment and growth of Fucus can also be inhibited by L. littorea on smoother surfaces or if snail densities are relatively high (Lubchenco 1983). These sorts of grazer habitat interactions have been demonstrated in a number of locations. For example, grazing by littorinids can inhibit algal canopies on mussel beds (Petraitis, 1987;Janke, 1990;Wilhelmsen and Reise, 1994). At different heights on the shore, or in different regions, the impact of grazing littorinids on community structure may be limited (e.g., Janke, 1990;Lindegarth et al., 2001). Observations of L. littorea suggest potential roles for mechanical or chemical properties of algae in influencing grazing behaviour (Watson and Norton, 1985a,b), although Barker and Chapman (1990) caution that it is not possible to extrapolate from single grazer/algae interactions to processes that structure the algal canopy on a shore.
Formal stock assessment and fisheries management has not been applied to L. littorea. Chapman (2006) considered the fishery to be completely unregulated in Scotland and a similar situation exists in Ireland (Cummins et al., 2002). Management measures in the 10 regional Inshore Fisheries and Conservation Authorities (IFCAs) in England vary, but some include a summer closed season (May-September) for winkle picking and a few authorities suggest size limits (no removal of snails small enough to pass through a 16 mm gauge).
The small economic and physical scales of the harvest, along with difficulties in defining the exploited stock and establishing a stock-recruit relationship help to explain the absence of fishery management models for L. littorea. What can be established are population densities on shores and shell size classes (e.g., Williams, 1964;Cummins et al., 2002). Size classes can be used as the basis for population models (e.g., Grady and Valiela, 2006). The current study therefore defines a size-based model for L. littorea and suggests harvest levels for shores based on the level of recruitment. A simple size-based model will produce a stable size structure. Such an outcome is not realistic for L. littorea, where size structure appears to change with density. Density-dependent processes are also known in other intertidal gastropods (Branch, 1975;Boaventura et al., 2003). As density dependence affects the responses of populations to harvesting (Bardos et al., 2006), the model for L. littorina estimated the potential for density to affect growth and survivorship. There are no specific restrictions placed on winkle picking within the lough and, as yet, no specific impacts of harvesting on features of conservation importance which would require intervention from the Northern Ireland Environment Agency. There is evidence of shellfish harvesting in both prehistory (McErlean et al., 2002) and more recent times. In the early twentieth century, 3 tonnes of winkles were reported to be harvested annually from Strangford for export to markets in Glasgow and Liverpool (Kelso and Service, 2000).
Population surveys were made over a total of 18 sites around the lough (Fig. 1). At each sampling event, 0.25 m 2 quadrats (n = 24) were haphazardly thrown, stratified into three contiguous bands down the shore: an upper band from the upper limit of Fucus spp. and L. littorea to the upper limit of Ascophyllum nodosum, a middle band where Fucus vesiculosus and A. nodosum are found and a lower shore band where Fucus serratus is found, bounded by the low water mark. All L. littorea individuals were removed from quadrats and length of the shell was measured. Larvae are thought to recruit to shores in the first half of the year (Moore, 1937;Williams, 1964;Fish, 1972), the smallest individuals in samples were therefore likely to have recruited to populations in the preceding year. The vertical stratification of sampling follows the methodology used by Cummins et al. (2002) in a survey of shores in Ireland. The horizontal extent of sampling for each site was approximately 30 m. All quadrats within a site were combined for population estimates as individuals have been inferred to migrate vertically as juveniles, seasonally, or in response to environmental stress (Smith and Newell, 1955;Lambert and Farley, 1968;Underwood, 1973;Gendron, 1977;Warner, 2001). Sites were surveyed in spring of 2004, 2005 and 2006. Not all sites were surveyed every year and at sites where the shore was large enough, two separate vertical transects were made. Observations on the ease of access to sites were quantified using the number of minutes required to move from a vehicle to the shore (Tab. 1). The shortest time (5 min) implies a convenient car parking space adjacent to the shore. A time of 10 min implies a longer walk, including cases where the shore is next to the road, but parking on the road is not practical. The longest times reflect extended walks across stretches of shore to reach a site, or that the access is across private land.
A comparative sample of L. littorea sizes at the point of sale was made from material bought at St. Georges market in Belfast during June 2005. Comparative data for mean shell lengths and densities are available for a survey of Irish shores made by Cummins et al. (2002) using the same methodology. Shell length and weight measurements were also made from a sample of L. littorea made at Silverstand in Galway.

Population characterization and modelling
The shell lengths of all field-sampled individuals from Strangford Lough (n = 11,408) were resolved into separate size classes for population modelling. The overall distribution of individuals was considered to consist of a small number of size classes where each size class is represented by a normal distribution of lengths. The unknown mean and standard deviation for each size class were estimated from a mixture model, which finds suitable parameters using an expectation maximisation algorithm. Different numbers of size classes were fitted using the normalmixEM function in R (Benaglia et al., 2009). The most parsimonious number of size classes was chosen after inspection of mixture model fits. We follow previous studies in assuming that size classes can be resolved as indicators of age in years (Moore, 1937;Williams, 1964).
Population densities projections were based on a sizestructured transition matrix model (Caswell, 2006) with density-dependent transitions from the first to second size class and from the second to the third size class.
Here L x,t is the density of L. littorea in age class x at time t (years), r is the annual recruitment of individuals from the plankton into the first age class, f(L x,t ) is a function describing the transition from size class x to size class x þ1. Given a value for r, the transition rate functions can be entered into a spreadsheet to predict densities in different size classes and the total density.
The recruitment rate is considered to be a constant. This reflects the open nature of populations, where the larval dispersal range is greater than the scale of individual sites (cf. Hyder et al., 2001). Density dependence was modelled using the Ricker curve to define the transition between size classes: Here the constants r and k have interpretations as intrinsic growth rate and carrying capacity. Theoretical treatments of resource competition can be used to derive the Ricker curve, providing potential mechanisms for density dependence (Geritz and Kisdi, 2004). Applying the Ricker curve on successive size classes implies that the most important density-dependent processes occur within a single size class (e.g., the transition from size class 2 to 3 is independent of densities in size classes 3 and 1). This simplification is consistent with the competitive effect demonstrated for Littorina close to 1.4 cm shell height (Petraitis, 2002) and the lack of a competitive displacement of other size classes when the densities of larger L. littorea were increased (Fenske, 1997). A higher density smaller size association, consistent with density-dependent processes, was found across the two sites investigated with Eschweiler et al. (2009). The higher density site also had evidence for higher rates of densitydependent mortality from parasites and crab predators (Eschweiler et al., 2009). Ricker curve fits were made to plots of the density in adjacent size classes at the same transect. This assumes that recruitment and mortality at each shore is not sufficiently variable from year to year to obscure the underlying transition rates.

Results
The distribution of L. littorea shell lengths was resolved into three classes by mixture modelling (Fig. 2). Transitions between size classes occurred at 1.1 and 2.06 cm, with mean shell lengths for each size class of 0.81, 1.56 and 2.48 cm. The average size of periwinkles at the fishmarket was 2.28 cm. The majority (69%) of winkles bought at the market could be considered as belonging to largest size class defined from field surveys from Strangford Lough.
Recruitment to shores was relatively predictable, with the shore-specific densities in the smallest size class tending to be correlated in separate (2004 and 2006) surveys (r s = 0.77, p < 0.05, n = 14). Therefore, where a shore had relatively high recruitment rank, it tended to retain this from year to year. The fraction of size class 1 (0.81 cm mean) L. littorea in the next size class declined with density (Fig. 3a). The Ricker curve fit to this transition parameters were r = 0.94 (SE. 0.146), k = 116.61 (SE 15.152), r 2 = 67%. The relationship between size classes 2 and 3 appeared to be noisier than the class 2 to 1 relationship. A Ricker curve emphasized a positive relationship between densities in adjacent size classes up to around 10 individuals m À2 . Parameters for the second Ricker curve were r = 1.19 (SE. 0.199), k = 13.70 (SE 1.419), r 2 = 19%. Once the transitions between year classes are defined, the total population density and size structure can be defined for any population, given a value for recruitment into year 1. As the rate of transition to larger size classes appears to decline with density, this implies that higher density populations will have proportionally fewer large individuals. This pattern was evident at 30 m transect scale and when pooling surveys by site (Fig. 4a). Predictions from equation (1) for the proportion of year 3 (> 2.06 cm) individuals were made by solving for different values of the recruitment rate r. These predictions agreed well with the observed values (observed-predicted regression F 1,16 = 75.85, p < 0.05, r 2 83%). In contrast, measures of shore accessibility to potential harvesters were not good predictors of total density (p > 0.05, r 2 5%) or the proportion of year 3 individuals ( Fig. 4b; p > 0.05, r 2 8%). Recoding the sites marked as private to  5 min access (in case landowners facilitate harvesters using private roads) did not change the results, and the strength of correlation between access time and proportion of larger winkles remained weak and not statistically significant.
A second prediction of the size-based model specified by the functions in Figure 3 is that mean size will decline with density. This decline in mean size was observed as a general trend in the Strangford survey and in data from the Cummins et al. (2002) survey (Fig. 5). A predicted mean-density relationship can be derived for populations described from equation (1). The predicted mean-density relationships for mean size on the basis of total density were realistic for both Strangford (observed-predicted regression F 1,49 = 56.34, p < 0.05, r 2 55%) and the more extensive Irish survey (F 1,52 = 7.12, p < 0.05, r 2 12%).
Larger periwinkles are more valuable at market. A simple harvest rule would therefore be to restrict sales to year 3 individuals, above 2.06 cm in length. This produces a peak in production of 13.8 year 3 individuals per m À2 at a year 1 density of five individuals m À2 . This would represent an annual harvest weight of 67 g m À2 , based on a mean weight of 4.82 g for year 3 individuals. The average weight for year 3 individuals given is based on an integration of the cohort density based on lengths (Fig. 2) and the length-weight relationship from Silverstrand (weight (g) = 0.349 Â length (cm) 2.854 , F 1,66 = 6233, p < 0.001, r 2 = 99%).

Discussion
The aggregation of harvests across different shores and a planktonic larval stage make it difficult to define a coherent L. littorea stock for management purposes and conventional fisheries modelling. The approach followed in this paper is to define a production rate per unit area, given a specified level of recruitment. The consistent ranking of year 1 densities across shores suggests that it is reasonable to describe shores by their level of recruitment. Size class definitions are similar to Williams (1964), who defined shell lengths at the end of the first year as 1.41 and 2.24 cm at the end of the second year (using length = height Â 1.66-0.04 derived from Silverstrand individuals to convert from shell heights in Williams 1964). The slightly lower size class boundaries in the current study (1.1 and 2.06 cm) may reflect spring surveys as opposed to Williams' estimates of size for the oldest individuals in a cohort. Due to the influence of density dependence, shores with the highest numbers of large individuals have intermediate recruitment.
Harvesting only the larger L. littorea is a process of passive management. This approach may be practically achievable by ensuring that wholesalers have appropriate grading equipment (Cummins et al., 2002). Larger L. littorea fetch a higher price, so market forces also incentivise focussing on larger individuals. If pickers are not paid for undersize winkles, hand picking will swiftly adapt to only picking the larger individuals. As L. littorea can become mature at 2 years old (at 1.1 cm shell height Williams 1964, equivalent to 1.83 cm length), harvesting larger individuals does not remove all the reproductive potential of a region. Furthermore, only accessible shores are harvested, areas that are not economic to visit could act as refugia for older reproductive individuals.
A reduction in the density and/or size of intertidal molluscs with human impact has often been observed (e.g., Keough et al., 1993;Roy et al., 2003). This pattern was not seen in the data from Strangford Lough. The lack of a relationship between a proxy of collection effort and both L. littorea density (Pearson's r = À0.1, p > 0.05) and the proportion of large (year 3) snails (Fig. 4), however, reflects patterns observed elsewhere. Tinlin-Mackenzie (2018) compared three shores of different collection effort in Northeast England, but neither L. littorea density nor median shell size clearly reflected collection effort. Unpublished reports cited by Tinlin-Mackenzie (2018) also had contradictory relationships between collection effort and size. Crossthwaite et al. (2012) reported smaller L. littorea at two high-harvest sites, but density was also higher at the same sites. As size covaries with density, the size differences observed by Crossthwaite Unfortunately, Crossthwaite et al. (2012) did not count L. littorea below 15 mm shell height, so their observed site differences cannot be interpreted using the size-structured model. In contrast, a reinterpretation is possible for the reserve effect reported on the shores of Gruinard (Johnson et al., 2008) as density and shell sizes are both available. Gruinard is a small island used for biological warfare tests during the Second World War. As a result of this, there were no visitors while the island was quarantined, and few visit now that the island has been returned to its original owners. A lack of access to harvesters was interpreted to be the cause of relatively high proportions of large L. littorea on Gruinard by Johnson et al. (2008). Mean shell lengths for two shores on Gruinard were 1.84 and 2.28 cm for densities of 121.7 and 4.3 individuals (all sizes) m À2 . Both these shell lengths exceed the predicted line in Figure 5. The pattern is opposite for two shores near Gruinard accessible to harvesters: mean lengths of 0.73 and 1.48 cm for densities of 50.7 and 35.8 L. littorea m À2 , respectively. The size-density relationships on and in the vicinity of Gruinard are therefore consistent with a reduction in the proportion of large individuals with harvester access. The absence of a pattern across Strangford Lough sites may reflect that access time is not a good proxy for harvest activity, or harvests may be too small or variable to have a detectable impact.
Having a basis for estimating the annual production of 3-year-old L. littorea has potential benefits for licencing initiatives. If the maximum harvest is 67 g m À2 , this is potentially a first sale value of €0.14 m À2 . A realistic shore segment of shore in Strangford Lough could be a harvestable area of 6000 m 2 , representing a potential annual harvest of €820. If a restricted number of harvesters are given access to a shore, a licence cost that gives an economic return per licence holder can be devised. An area-based harvest is also applicable to a number of management frameworks, including Territorial Use Rights in Fisheries (TURFs), where small-scale fishers collectively manage the stock (Leiva and Castilla, 2002). As harvesters are thought to target larger individuals, a weight or bag limit could be devised based on the number of harvesters and the predicted weight of snails in a region. Even in the absence of a licencing or management framework, a basis for predicting the likely harvest on different shores can be useful in quantifying ecosystem services and understanding trade-offs in marine spatial planning. Given a minimum of assumptions, the likely harvest for any particular area can be estimated using the model and compared to catch estimates for a preliminary evaluation of the sustainability of harvests.
The model could also be used to evaluate different management approaches. For example, thinning of year 2 individuals on some shores could increase the production of older snails. This would be a good test of the model's mechanism, but given the relatively low value of the harvest, additional effort in the field is unlikely to be economically sustainable. Other management options are less easy to evaluate, regardless of a model. Having a closed season does not automatically result in a reduction of effort. If a seasonal closure is to conserve stocks, a reduction in harvest would need to be demonstrated before any estimates of impact could be made. Similarly, closed areas may increase the potential supply of recruits to shores, but evaluation of such an approach would need a combination of hydrodynamic modelling, surveys of recruitment and population genetics. One prediction of the population model developed in this paper is that, even if protected areas increase recruitment to a shore, the population of year 3 individuals may not increase due to the observed density dependence. A similar conclusion is in the marine protected area (MPA) literature: that density dependence can explain why some reserves do not result in increased fishery yields (Gårdmark et al., 2006).
The advantage of a model is that harvests can be predicted even if all year 3 individuals have been removed at a site. The Cummins et al. (2002) data were not used in deriving parameters, indicating that the model may have some generality. The actual level of harvesting across the Strangford and Cummins et al. (2002) survey sites was not quantified. This is likely to have added to the uncertainty in parameter estimation. Potentially, some of the Cummins et al. (2002) survey sites may have a lower level of harvesting, reflected in apparently higher mean sizes when compared to the Strangford survey. A potential issue with the model is ignoring ages above 3 years. This may not be significant as such individuals seem rare in the surveyed populations. Time series of survey and catch from individual shores would help to evaluate any biases in the single survey-based model presented in this study. A further consideration is that parameter values may vary with factors like wave exposure or climate. It would be possible to tune the model for additional environmental covariates using more extensive data. Similarly, the transition functions of the model could be modified by climate change, for example as growing and recruitment seasons alter. Simkanin et al. (2005) considered L. littorea to be a northern species and recorded a decline in abundance between surveys of the Irish coast made 45 years apart. On the basis of fisheries figures, Simkanin et al. (2005) interpreted the L. Littorina decline to fishing, not climate. Such an interpretation may not be justified given the low confidence in catch figures for L. littorea and the market forces that may also affect harvest activity. Transitions in the model show nonlinearity, but also the initial slopes are above 1, indicating that at low densities, more year 2 individuals are found than year 1 individuals (and initially year 3 density > year 2). This can be interpreted as immigration at low densities and/or not all individuals being found in surveys. Undetected recruits seem likely considering the long recruitment season and the small size of the initial benthic stages. It may be possible to find ways to sample more completely, but the surveys were carried out with one shore per tide and it would probably not be practical for management to develop a more time-consuming approach.
The size-structured model for L. littorea is not an "integrated" model as it does not use multiple data types (Punt et al., 2013). In particular, there is no clear link to fisheries data, reflecting the aggregated and unreliable landings data. A key element of the approach is that density-dependent transition rates were identified. Density-dependent transition rates are rare in size-structured fishery models, but they appear key to understanding the structure of L. littorea populations. Density dependence provides an explanation for observed declines in mean shell size with density. Considering the interrelationship of size and density is likely to provide clearer evaluations of harvesting effects than approaches that have tended to ignore this covariation. The L. littorea model can be used to estimate the potential catch from harvested shores, adding in valuation of resources and guiding potential catch limits. As the approach uses size and count data, it is relatively simple to redefine the parameter values for other locations using additional surveys. The integration of size class and site data in the current study may provide a framework for studying other data poor spatially dispersed stocks.