Spatial distribution patterns and population structure of the sea urchin Paracentrotus lividus ( Echinodermata : Echinoidea ) , in the coastal fishery of western Sardinia : a geostatistical analysis

The identification of appropriate spatial distribution patterns for the observation, analysis and management of stocks with a persistent spatial structure, such as sea urchins, is a key issue in fish ecology and fisheries research. This paper describes the development and application of a geostatistical approach for determining the spatial distribution and resilience of the population of the sea urchin Paracentrotus lividus in a fishing ground of western Sardinia (western Mediterranean). A framework combining field data collection, experimental modelling and mapping was used to identify the best-fit semivariogram, taking pre-fishing and post-fishing times into consideration. Variographic analyses indicate autocorrelation of density at small distances, while the isotropic Gaussian and spherical models are suitable for describing the spatial structure of sea urchin populations. The point kriging technique highlights a generally patchy population distribution that tends to disappear during the fishing season. Kriging maps are also useful for calculating predictable stock abundances, and thus mortality rates, by class diameters within six months of fishing. We conclude that the framework proposed is adequate for biomass estimation and assessment of sea urchin resources. This framework can therefore be regarded as a useful tool for encouraging a science-based management of this fishery.


INTRODUCTION
Over the past decade there has been increasing interest in modelling and measuring spatial patterns (e.g.gradients and patches) in biotic variables as a means of understanding the mechanisms that control critical aspects of the ecology of species, such as spatial distribution (Legendre and Legendre 1998).
With the advancement of computer science, geostatistics has become a powerful tool for estimating the spatial distribution of marine populations (Conan 1985, Maynou 1998), for predicting stock abundances (Petitgas 1993, 2001, Megrey and Moksness 2009) and for assessing marine reserve benefits (Stelzenmüller et al. 2007).Indeed, the identification of appropriate spatial scales for the observation, analysis and management of stocks with a persistent spatial structure is a key issue in fish ecology and fisheries research (Orensanz et al. 2006, Ciannelli et al. 2008).
Owing to their low mobility, numerous benthic commercial species can be considered suitable for geostatistical applications (Jensen andMiller 2005, Adams et al. 2010), and the purple sea urchin (Paracentrotus lividus) is an ideal species for a case study.P. lividus is common throughout rocky intertidal and shallow subtidal zones of the Mediterranean Sea and northeastern Atlantic Ocean, where it is generally associated with erect macroalgae.Its spatial distribution can vary on both small and large scales in relation to the interaction of abiotic and biotic factors (Boudouresque and Verlaque 2001).For example, the variability of water temperature and solar radiation can partially explain intraspecific variation in covering behaviour type in P. lividus (Crook 2003).The heterogeneity of the substratum plays a key role in providing P. lividus with shelter, thus influencing the structuring of populations, where predation pressure (which includes human harvesting) is particularly high (Bonaviri et al. 2005, Hereu et al. 2005).Among human-related impacts, site accessibility during harvesting by diver fishermen significantly affects the structuring of sea urchin populations in a fishing ground in northern Sardinia (Ceccherelli et al. 2011).Such studies, like other significant ones on sea urchin predation (Sala andZabala 1996, Guidetti et al. 2004), recruitment (Tomas et al. 2004), migration (Palacín et al. 1997, Crook et al. 2000), competition (Guidetti 2004) and harvesting (Pais et al. 2011), employ conventional approaches that assume spatial independence of a measured variable (specifically abundance indices), i.e. values at one location are independent of values at neighbouring locations.Although conventional approaches are equally valid, they involve some limitations in the usefulness of the ecological data gathered for biomass estimates, in terms of spatial scales and the setting of their confidence limits (Addis et al. 2009).In this respect, geostatistical techniques are more powerful tools for estimating the spatial distribution of marine benthic communities than conventional statistical methods because they explic-itly consider spatial correlation between observations (Warren 1998, Rueda 2001).
According to Boudouresque and Verlaque (2001), the spatial domain of P. lividus populations ranges from fishing areas where the stock is sufficiently abundant to support a commercial fishery (an isolated bay of a few square kilometres) to small-scale aggregations within a bed or "patches", measuring tens to hundreds of square metres, where ecological experiments are usually carried out.However, there have been no attempts to describe the spatial structure of P. lividus populations by geostatistics, for instance by estimating semivariograms and their descriptors (nugget, sill and range), which are useful for evaluating the extent of spatial correlation in the data.Such applications are the basis for spatial perception of sea urchin stocks and thus for the success of fisheries management (Chen andHunter 2003, Grabowski et al. 2005).
P. lividus is the main echinoid exploited in Europe (FAO 2011), but the information we have on the current status of populations in the coastal areas of the Mediterranean is scant (Andrew et al. 2002).Trends in relative abundance or stock assessments have never been estimated, though the decline of landings in a few fisheries indicate that populations have been severely depleted (Andrew et al. 2002).Major concerns regard the lack of time series data in terms of both commercial landings and fishery-independent surveys, which are indispensable for the assessment of sea urchin stocks by means of catch-per-unit-of-effort (CPUE) (Perry et al. 2002, Chen andHunter 2003).
We underline that sea urchin fisheries for P. lividus are in need of a precautionary approach "sensu FAO" (1996) in order to avoid the risk of stock collapse, as has occurred in some fisheries of northern Europe, where studies on the impacts of sea urchin harvesting were neglected (Sloan 1985, Byrne 1990).Our case study refers to the sea urchin dive fishery of Sardinia (southern Italy), where fishing for P. lividus has a significant social-economic impact but the management has been largely unsuccessful at conserving the stock and ensuring a sustainable fishery (Pais et al. 2011, Cau et al. 2007).
Since the elucidation of spatial distribution patterns is essential for abundance estimates of sea urchins stocks, the objectives of this study were to a) determine, model and map the spatial structure of the P. lividus population in a fishing ground in western Sardinia; b) assess the spatial patterns by size considering the pre-and post-effects caused by fishing harvesting; and c) evaluate the predictable number of specimens to detect abundance fluctuations due to harvesting and assess the total mortality rate.

Study area and field samplings
The study area is located in central-western Sardinia (Fig. 1).Cape Pecora is a shallow open bay where the geomorphology of the sea bottom is characterized by pebbly, metric and decimetric blocks that are highly eroded and rich in caves and rock shelters, providing a suitable habitat for certain benthic organisms, particularly sea urchins.The biocenosis includes upper subtidal algae exposed to high wave action (UIW) and Posidonia oceanica meadows (Pérès and Picard 1964).The most representative algae are in the genus Cystoseira, with a prevalence of Cystoseira stricta var.amantacea.Other important algal genera are Laurencia, Dictyota, Dictyopteris, Codium, Stypocaulon, Padina, Acetabularia, Halimeda and Amphiroa.
Surveys were conducted in an area encompassed by a 1.5-km stretch of shoreline to a depth of 10 m (with a total surface area of 0.2 km 2 ).We superimposed a regular grid, subdivided into 30×30 m cells, and selected 90 of the grid's 270 cells, representing one-third of the whole area.Stations were randomly selected as starting points from which underwater counts of sea urchins were carried out within three random replicate quadrates of 1 m 2 .Each station was geo-referenced (latitude-longitude) by GPS using Universal Transverse Mercator projection (UTM).Size data were obtained for each station by measuring all individuals larger than 1 cm in diameter using a Vernier calliper (mm).Data on diameter size were grouped into three size classes: 10-29.9mm (Juvenile), 30-49.9mm (Medium) and ≥50 mm (Adult), which corresponds to the minimum size for commercial fishing.
Experimental surveys were conducted in October 2010, prior to the beginning of the fishing season (pre-fishing) and in May 2011 at the end of the fishing season (post-fishing) for all of the 180 stations investigated.

Statistical analysis
Mean density (±SE) for three size classes was calculated; results were plotted on histograms and classed post-maps, to check for errors in the raw data and to verify whether geostatistical analysis could be applied.A preliminary Z-test was performed to detect differences between mean pre-fishing and post-fishing density.
Density was successively defined as degrees of autocorrelation between measured data points for each size class diameter.This was obtained by a non-directional experimental semivariogram γ(h) computing the variance of a population, while taking the spatial position of the sampled stations into account, making use of the equation (Matheron 1965 where Z(x i ) represents the density of sea urchins at sampled station x i , Z(x i +h) is a variable value separated from x i by a distance h (measured in metres), and N(h) is the number of pairs of observations separated by h.To avoid decomposing semivariograms at large lag intervals, the default active lag distance was set close to 70% of the maximum lag distance.We undertook the following estimation for each experimental semivariogram: the nugget effect (C 0 ) is attributable to measurement error, micro-scale variability or small-scale spatial structure; the sill (C+C 0 ) can be defined as the maximum variability point beyond which the semivariance values become asymptotic; the range (A 0 ) represents the distance within which the data remain autocorrelated (Maynou 1998).The model that best explained the spatial structure of each case was selected on the basis of values for the reduced sum of squares (RSS) and the coefficient of determination (r 2 ) (Cressie 1991).
Semivariogram parameters of the selected model for each size class at each time were employed using the spatial estimation technique known as "point kriging".This enabled us to create two-dimensional density maps.The kriging estimate of Z(x) at each node was obtained by a linear combination of the samples, each weighted by a factor (λ i ), which depends on the combination of the relative position of the sampling points, the theoretical semivariogram, and the Z(x i ) values at the sampling points (Matheron 1965).The estimated density values Z(x) were given by The validity of the models in the variographic analysis and kriging interpolations was evaluated using jack-knife cross-validation, performed by sequentially deleting one datum and using the remaining data to predict the deleted density value; the selected semivariogram model and kriging parameters were applied to this end (Maravelias et al. 1996).The observed (O) and estimated (E) densities were plotted and fitted to a linear regression O=α+βE; the significance of α and β was tested (t test) under the null hypotheses α=0 and β=1 (P=0.05)(Power 1993).
The predictable number of specimens was calculated by scaling the surface of kriging maps with mean densities for each countering layer, including the confidence limits (mean±SE).Predictable number of specimens of pre-fishing (October) and post-fishing (May) by size was used to estimate the total mortality rate (Z=M+F) (Ricker 1975) by Z = -ln N t /N 0 where N t is the estimated number of sea urchins in the post-fishing period and N 0 is the number in the prefishing period.Since Juveniles and Medium individuals should only be affected by removal not associated with fishing, the total mortality rate (Z) for these classes corresponds to the natural mortality (M).
Calculations of semivariograms and kriging maps were carried out using Gs+ ver.7 (Gamma Design Software, LLC) and Surfer8 (Golden Software, Inc.) geostatistics software.
The calculation of experimental semivariograms revealed that densities of the three size classes of P. lividus were spatially structured (Fig. 3).Isotropic Gaussian and spherical models provided the lowest RSS and the highest r 2 of all analysed models, successfully explaining the spatial population structure of the three size classes analysed from the two periods (Table 1).In all cases, there was no significant discontinuity at origin (C 0 ≤19%), thus indicating that the sampling spa-  tial resolution used was appropriate (Fig. 3).In the prefishing phase the semivariogram of Medium provided the highest range value (A 0 =132.00m), while that associated with Adult was the lowest (A 0 =33.20 m); a similar result occurred in the post-fishing period but in this case the size class of Juvenile showed an increase, with the distance of spatial influence varying from 48.32 to 100.8 m (Table 1).All the size classes showed the lowest values of sill in the post-fishing period (Table 1).The spatially structured density component [C/ (C 0 +C)] varied between 81.70% and 99.90%, indicating that the sea urchin population that we studied had a well-defined spatial structure (Table 1).
Cross-validation analysis supported the appropriateness of the Gaussian and spherical models and kriging predictions; in all cases the values tested by the t test revealed no significant differences (P>0.05).
The densities estimated by kriging allow for visualization of spatial distribution by size over time (Fig. 4).During pre-fishing, the Juvenile density remained within the range of 0-2 ind.m -2 throughout the study area.A few sporadic "density hot-spots" (with densities up to 5 ind.m -2 ) were localized along the coastline.Individuals belonging to the Medium class showed a patchy distribution, with densities generally ranging between 1 and 3 ind.m -2 , but with some hot-spots where densities higher than 4 ind.m -2 were recorded.The density of Adults was consistently low (<2 ind.m -2 ), with only hot-spot density reaching levels as high as 3 ind.m -2 (Fig. 4).
Densities in the post-fishing phase were lower than those in the pre-fishing phase.In details, Juveniles were missing from the northern area and occurred in low numbers, never exceeding 4 ind.m -2 in the southern areas, with only two isolated higher-density hotspots appearing in the post-fishing map.Densities in the Medium class showed a similar pattern to that of the Juvenile class, but in the latter case the hot-spots had almost disappeared, with only one case (of 5 ind.m -2 ) occurring in a large area.The density in the Adult class ranged from 0 to 3 ind/m 2 , with a patchy distribution during the periods studied.The number of patches was also noted to be lower than that measured during the pre-fishing period (Fig. 4).
The predicted numbers (mean±SE) of Juveniles in the pre-and post-fishing periods were 535274±37016 and 408331±31879, respectively, indicating a total mortality rate Z=0.27.The estimated numbers for the Medium class in the pre-and post-fishing periods indicated values of 1119652±68854 and 907796±59 005, respectively, with Z=0.21.The ranges of estimated number of harvestable sea urchins (Adult) in the pre-and post-fishing periods were 19068±1580 and 11701±951, respectively, indicating a total mortality rate Z=0.49.

DISCUSSION
Through the use of geostatistics we have described the patterns of spatial distribution of the purple sea urchin in an area of the western coast of Sardinia.In this region the commercial harvesting of sea urchin is commonly practiced by the dive fishery.The case study thus reflects a true condition of an exploited stock.The geostatistics application and the kriging maps explain the spatial characteristics of the resource on both spatial and time scales.Density maps were useful for predicting the likely stock biomass by area, which is important key information for the development of spatial management-based quotas.
Results of variographic analyses show autocorrelation of density at small distances.The isotropic Model parameters are listed in Table 1.
Gaussian and spherical models describe the spatial structure of sea urchins in the area, confirming a previous study conducted in Sardinia (Addis et al. 2009).
We observed a general reduction of semivariance and changes in the spatial structure of sea urchins from the pre-fishing to the post-fishing, which is a consequence of the decrease in sea urchin densities owing to fishing over a six-month period.
The most important changes were related to spatial patterns for the Juvenile class.Indeed, the kriging map in the pre-fishing phase shows a clear hot-spot pattern with a regular spread of contour levels (Fig. 3).This result is particularly interesting because it highlights the existence of nursery sites as a result of the good fit of settlers for the seabed.In the study area the geomorphology of the bottom is characterized by caves and rocky shelters, which provide protection for small sea urchins that are "post-settlers" from predators and destructive wave action, the main perturbation affecting populations (Sala et al. 1998, Hereu et al. 2004).In the scientific literature, attempting to explain the settlement process of in sea urchins has raised conflicting opinions.Some authors have pointed out that there is an absence of significant autocorrelations between sampling sites located even tens of metres apart, demonstrating a high level of settlement heterogeneity over very small patch sizes (tens of metres) (Hereu et al. 2004).On the other hand, other authors claim that settlement of echinoids only shows spatial variability on a scale of thousands of metres rather than on smaller scales (Keesing et al. 1993).Considering that the home-range of sea urchins is generally low, varying from 50 to 300 cm (Hereu 2005), we can assume that the Juveniles were "new settlers".Our results concurred with those of Hereu et al. (2004), who found high spatial variability over small distances and very small patch sizes.Analysis of the post-fishing phase for Juveniles suggested a change in the spatial pattern.Considering that this class would not be affected by fishing harvesting but only by natural mortality, specimens from this class would be expected to be distributed over the whole area.Nevertheless, our results indicated that these specimens disappeared from many areas and were confined to the southern area.We propose two assumptions to explain this observation.First, that these specimens had shifted into the Medium class as a consequence of growth during the six-month experimental period, in which case we would expect an increase in diameter lengths of 4 mm based on the growth parameters proposed by Allain (1978).Our second statement explains the pattern identified as a consequence of mass-mortality caused by strong events of Mistral wind (NW) during the studied period.This process, which is common in shallow, exposed habitats (Turon et al. 1995), might affect the Juvenile class and thus the population size structure and density over time.
Sea urchins belonging to the Medium class indicated differences in spatial distribution between pre-fishing and post-fishing.As noted for the previous class, the semivariogram described the autocorrelation of data changes from the spherical to the Gaussian model.The main changes encompass a decrease in the estimated number of sea urchins (from 1119652±68854 to 907796±59005) and the occurrence of only one density hot-spot.Besides natural mortality and the shift into the next size class, this group could also be subjected to fishing mortality because specimens are illegally harvested by commercial fishing operators (Cau et al. 2007).
The spatial structure for large specimens (Adult) did not show significant changes over time.A decline in semivariance of this group was due to a decrease in the population densities, based on densities measured during the pre-fishing and post-fishing periods.The patchy distribution identified was relatively stable over time but there was a substantial decrease in the estimated number of sea urchins (from 19068±1580 to 11701±951).In this case the total mortality for the Adult group was quite high (0.49): double that of the Medium class.At this size sea urchins are affected by natural and fishing mortality (Z=M+F).Considering a steady condition, we expected to identify a total mor-tality similar to that of the previous classes (0.27 and 0.21, respectively).Since we observed a total mortality rate of 0.49, it is very likely that fishing mortality is approximately 0.2.
Reduction in density, size structure and biomass of P. lividus as a consequence of fishing has been reported in the literature (Pais et al. 2011, Guidetti et al. 2004) but the method proposed in the present paper provides a useful technique that permits the quantification of these reductions, highlighting where and how they occurred.Furthermore, kriging maps are easier to read than other statistic outputs, and can therefore become a useful tool for the main stakeholders tasked with managing the sea urchin dive fishery.
Spatially explicit management strategies are a key requisite for maintaining persistent stocks, such as sea urchins, that are dominated by low mobility and fragmentary and restricted fishing grounds.According to scientific surveys conducted on the whole of the Sardinian coast in 2006-2007 (Cau et al. 2007), the domain of sea urchin population is on a scale of a few kilometres (mesoscale).This spatial scale corresponds to the domain of the sea urchin dive fisheries that harvest the resource in fixed fishing grounds.The assessment of sea urchin stocks around Sardinia using the model we tested is thus a pragmatic tool for predicting abundances per fishing ground.In terms of developing a scientific-based management of the sea urchin fishery (Perry et al. 2002, Chen andHunter 2003), we stress the development of a systematic scientific programme aimed at acquiring substantial data that should be interpreted together with fishery-dependent statistics.

Fig. 1 .
Fig. 1. -Study area at Capo Pecora bay (W Sardinia, Mediterranean Sea) showing the fixed grid of 270 cells and the 90 sampling stations (•) randomly selected during the pre-fishing survey.

Fig. 2 .
Fig. 2. -Percentage size distribution of sea urchins in pre-fishing and post-fishing.

Fig. 3 .
Fig. 3. -Experimental semivariograms estimated with spherical and Gaussian models for each size class during pre-fishing and post-fishing.Model parameters are listed in Table1.

Fig. 4 .
Fig. 4. -Kriging maps for P. lividus population at Capo Pecora bay based on geostatistical interpolation at pre-fishing and post-fishing.Mean densities of Juvenile, Medium and Adult urchins are shown.