The Pedro Bank is a shallow offshore bank located in the west-central Caribbean approximately 80 km south of the main island of Jamaica (Fig. 1). The bank has a total area of 8040 km² with maximum dimensions of 180 km west to east and 79 km north to south down to the 30 m depth contour (Morris 2016). The average depth is 24 m and the bottom consists predominantly of a coralline sand substrate which makes up approximately two-thirds of the bank. Infused throughout this primary substrate are patch coral reefs, coral pavement (weathered coral reefs), various species of macroalgae and gorgonians (Baldwin 2015). This habitat complex accommodates a marine ecosystem that supports important fisheries for various finfish, spiny lobster (Panulirus argus) and queen conch (Baldwin 2015). The central and southern areas of the bank are known to contain the highest densities of both juvenile and adult conch (Morris 2016).

Distribution data were compiled from underwater visual surveys conducted by the National Fisheries Authority of Jamaica (NFA) for the years 2007, 2015 and 2018. The NFA has conducted such surveys at three to five-year intervals since 1994, but only after 2007 did sampling target the entire extent of the bank. The years 2007, 2015 and 2018 were chosen to contrast sample design, sample sizes and spatial population structures such as the size and location of high-density patches. The 2007 survey consisted of 60 sites clustered in the centre of the study area, while the 2015 and 2018 surveys included 80 and 81 sites, respectively, systematically distributed down to the 30-metre contour (Fig. 2). Sample sites were defined as one to four belt transects with an area of 220 to 2250 m² spaced one metre from the end of one transect to the start of the following one. The number and area of transects were dependent on ocean conditions and diver safety considerations at the time of sampling. Data recorded within each transect included geographic location, depth, transect area, habitat type and numbers of adult and juvenile conch. Juveniles were defined as individuals of shell length lower than 20 cm and without a folded shell aperture or lip, while adults were defined as having a shell length of 20 cm or greater and having a thickened or folded lip (Appeldoorn 1988). In recognition of the different data assumptions of the spatial models, calculated densities were used for the geostatistical models OK and KED, whereas count data were used for the other models.

We collated broad-scale spatial environmental layers for geographic location (latitude and longitude), habitat type and depth (bathymetry) of the bank to be used as covariates in the interpolation models. Geographic location was selected as a proxy for interactions between queen conch location and ecological processes (Gutierrez et al. 2008). Juveniles and adults aggregate in specific habitats for feeding, reproduction and predator avoidance, thus displaying heterogeneous distribution with respect to habitat (Stoner and Ray-Culp 2000). The habitat layer used was a generalized substrate categorization created at 250 m² resolution and adjusted to 100 m² from remote sensing and field measurements as part of the Pedro Bank Marine Spatial Planning Project (Baldwin 2015). Because the queen conch also displays stratification with depth (Morris 2016), we extracted depth readings down to the 30 m contour, coinciding with the maximum depth of the surveys, from the General Bathymetric Chart of the Oceans Compilation Group grid bathymetric database (GEBCO 2020).

We created a 100 m² gridded interpolation surface down to the 30 m depth contour consisting of 606000 pixels for a total area of 6060 km². Because the sampling area for 2007 was smaller, the grid was adjusted to 204000 pixels or approximately 2040 km² to avoid extrapolation for that dataset. Each environmental layer was converted to the same 100 m² resolution so that the interpolation outputs were comparable. For the interpolation process, each pixel was considered a discrete unit, so the interpolated value in each pixel was multiplied by the ratio of the total area of one pixel (1 km²) and the mean site area for each survey year to compensate for the difference between the unit used for the interpolation surface (pixel) and the area of the sampled sites (1 m²). Total biomass within a pixel was calculated by transforming the estimated abundances using queen conch weight conversion factors corresponding to the 50% clean processing level (shell, viscera and opercula removed) for adults (Aspra et al. 2009) and 36.69 g per individual for juveniles (Appeldoorn and Rodriguez 1994).

We calculated two design-based (Cochran 1977) biomass estimations for comparison with the spatial interpolation methods: the first was a non-stratified method and the second one was stratified. The non-stratified method consisted of an estimate based on the product of the mean abundance and the total area and the weight conversion factors for adults and juveniles. Design method estimations can be significantly improved by stratification (Chang et al. 2017), so we introduced stratification based on depth (0-10, 11-20 and 21-30 m). The second design biomass was then calculated similar to the first method but for each stratum and then summed to obtain total biomass.

Ordinary kriging and KED were applied to the datasets using functions and procedures in the R gstat package (Pebesma 2004). Geostatistical models use the variogram to estimate the spatial autocorrelation and covariance structures of the sample data, which are then used to determine kriging weights (Webster and Oliver 2007).

where γ h is the semivariance at the distance interval h , Z x i is the number of conches at sample site x i and N ( h ) is the number of observation pairs separated by distance h .

Ordinary kriging (OK) assumes a stationary local mean around each sample and uses an unbiased linear combination of the surrounding weighted values produced by the variogram to estimate values at unsampled sites (Webster and Oliver 2007).

where Z O K is the estimated density of conches at the unsampled point of interest, x 0 , Z is the observed value at site, x 0 , λ i O K is the weight and is the value of samples within the search neighbourhood used for the estimation. We defined the local neighbourhood by designating a maximum distance over which semivariance values are calculated and the minimum and maximum number of observation points to be used in the kriging prediction. Defining a local neighbourhood in this way reduces bias from outliers and extreme values of highly skewed data (Surette et al. 2007).

KED assumes that the local trend or drift is a linear function of a secondary variable (Rivoirard et al. 2000), which in this case was the depth covariate.

where f x is the drift function of the depth variable and a and b are linear coefficients. The KED estimation Z K E D ( x 0 ) is achieved by imposing two conditions: one where ∑ i = 1 n λ i K E D = 1 as with OK, and one where ∑ i = 1 n λ i K E D f ( x i ) = f ( x 0 ) is the secondary variable at each sample site x i and f ( x 0 ) is the globally known secondary variable at points x 0 .

A spatial version of GAM (Hastie and Tibshirani 1990) assuming a negative binomial probability distribution (NBGAM) was applied using the R mgcv package (Wood 2011), while the mixed model GAM or GAMM was applied using the R nlme package (Pinheiro et al. 2021) with procedures from Zuur et al. (2009).

where Y i is the number of conches at site i , is the intercept parameter, s x i 1 … p are smoothing functions of the covariates and ε i is the error term for the unexplained information. The starting model for each dataset was specified to include all available covariates (geographic location, depth and habitat). The best set of covariates was selected through a backward selection process using the Akaike information criterion (AIC) .

The GAMMs included optimal variance and correlation structures to model the random error term , heterogeneity of the variance and residual correlation in the data. Optimal variance and correlation structures were selected for the GAMM by minimizing the AIC among competing structures (Zuur et al. 2009).

Two hybrid regression kriging (RK) models (Hengl et al. 2007), GAM+OK and GAMM+OK, were used. The RK prediction Z R K at the target location x 0 was reached by summing the trend estimated by the regression procedure m ( x 0 ) and the kriged residuals using OK

The spatial negative binominal zero-inflated model (ZINB) directly models the two sources of excessive zeros which can cause bias in spatially sampled data: sampling error and ecological processes (Martin et al. 2005). The procedure included modelling the non-zero part of the data assumed to be negative binomially distributed and the probability processes that generate zeros. The selection of covariates for both processes was done using backward selection among the set of explanatory variables (geographic location, depth and habitat) by minimizing the AIC. The analysis was implemented using the zeroinfl function of the R pscl package (Zeileis et al. 2008).

Statistical evaluation and comparison of spatial model performance were done using ten-fold cross-validation. Two diagnostic statistics were calculated, the mean absolute error (MAE), a measure of prediction bias, and the root mean squared error (RMSE), which reflects prediction precision or prediction error (Willmott 1982). A smaller MAE or RMSE value represents a better performing model.

where n is the number of samples, p i is the predicted value at a location and o i is the observed value.

The performance of the spatial models reflected their underlying interpolation assumptions and the structural characteristics of the survey datasets. There was no clear correspondence between bias measured by the MAE and the prediction precision based on the RMSE. In other words, a model returning the lowest MAE did not necessarily return the lowest RMSE (Table 2). Only for the 2015 adult conch estimations did the model with the lowest MAE and RMSE relate to the same model. The best-performing models in 2007 and 2018 were generally those that included GAM or GAMM in full or in part through regression kriging (NBGAM, NBGAM+OK, GAMM and GAMM+OK). For 2015 the best models also included OK and ZINB. Specifically for the juvenile data, the best models were ZINB (MAE=0.13) and GAMM+OK (RMSE=14.46) in 2007, GAMM (MAE=0.06) and ZINB (RMSE=15.17) in 2015, and ZINB (MAE=0.27) and GAMM (RMSE=23.50) in 2018. For the adult data, the best models were NBGAM (MAE=0.57) and GAMM (RMSE=20.00) in 2007, OK (MAE=0.23, RMSE=33.60) in 2015, and NBGAM (MAE=0.15) and GAMM (RMSE=22.59) in 2018.

Adult biomass estimated from spatial models varied approximately in a range equal to the minimum estimate for each year, with a CoV of 0.40 in 2007, 0.15 in 2015 and 0.24 in 2018. For the juvenile biomass, the CoV was 0.60 in 2007, and 0.18 in 2015 and 2018. This indicates that there were no great differences between the results of the models in terms of biomass estimation.

The design-based estimations tended to be close to the biomass estimated by spatial models (Tables 2 and 3), but a major limitation is that the simplicity of this biomass estimation is not enough to develop ecosystem-based management (Garcia et al. 2003). For adults, the percentage differences between the design methods and spatial interpolation biomass estimations ranged from -5% to -14% and +7% to +43% for the stratified and non-stratified methods, respectively. For juveniles, the differences ranged from +33% to -43% and +6% to +51%, respectively.

The spatial distribution maps of the model predictions (Figs 4 and 5) highlighted the concentration of adult and juvenile conch biomass along the southern edges of the bank, as has been reported by Morris (2016). This area is characteristically shallower and contains a more complex habitat structure, including coral reefs, seagrass, gorgonians and macroalgae, as opposed to predominantly sand and algae in the other areas. However, each model differed in its degree of local versus global smoothing and identification of the local and broader-scale spatial structures. For instance, the OK, KED and RK models better identified local high-density patches while the NBGAM, GAMM, and ZINB models identified larger-scale distributional trends largely reflecting the environmental covariates. High-density patches of juveniles were generally fewer and smaller in terms of abundance in comparison with adults. They were also farther away from the southern edges of the bank than the adults, suggesting a slight distributional difference. The spatial distribution maps therefore validated the performance of the interpolation model by identifying conch aggregations related to the environmental covariates.

Modelling the spatial component of exploited sedentary stocks may be critical to precise biomass estimations, maintaining reproductively viable population structures and avoiding the pitfalls of recruitment overfishing. In this study, the models that assumed a priori probability distribution or used covariates to model the spatial structure (NBGAM, NBGAM+OK, GAMM, GAMM+OK and ZINB) generally performed better than those that did not, such as OK. This contrast clearly shows the importance of closely considering the spatial component of the abundance index being estimated related to the species biology and exploitation. Queen conch, like many sedentary species, aggregate in specific habitats to carry out biological activities and are subject to high levels of exploitation which constantly alter their population structure (Stoner and Appeldoorn 2021). While the weight of evidence has long favoured model-based approaches over design-based ones, the set of models is often generic in their calibration or is otherwise limited by the availability of covariates (Li and Heap 2011, Chang et al. 2017). In our case, careful knowledge of the species biology, level of exploitation and the study area were required to select and specify the models that were used.

From our results, the design methods returned biomass estimations that were sometimes very close to the best interpolation models (Table 3). Despite this, the presence of unit stock assumptions and failure to account for the population spatial structure renders design methods inadequate for sedentary species (Anderson and Seijo 2010), so they should not be applied to queen conch abundance index estimation. In other cases, the design methods returned estimations with large differences from the best interpolation models. From an ecosystem management standpoint, the risk of estimation errors is particularly problematic for density-dependent sedentary species such as the queen conch, in which allee effects can be enhanced by fishing leading to negative population growth and stock collapse (Stoner and Ray-Culp 2000). Chang et al. (2017) suggested that predictions from design methods may be comparable to those from statistical models if the study area is accurately stratified. Accurate stratification can, however, be difficult or subjective and may contain large bias due to the unmodelled spatial autocorrelation from the spatial structure (Segurado et al. 2006). In our case, we stratified according to ocean depth, which generally improved on the non-stratified estimate in terms of the relative difference from the best interpolation model. Other stratification criteria could have also been considered, such as habitat type or age class, but these were either impractical to apply or the data were not available.

In this study, the best models were generally the GAMM and the ZINB based on the two evaluation statistics, MAE and RMSE. It was, however, clear that model performances were strongly influenced by the structure of each dataset, including its sampling design, sample size, coefficient of variation and the quality of the available covariates used. As a result, the best model was not consistent for either juveniles or adults across the three surveys. This is consistent with theory and studies that suggest that spatial interpolation model performance is often data- and case-specific (Isaaks and Srivastava 1989, Rufino et al. 2021) and is reflected in the level of bias and prediction accuracy (Li and Heap 2011). For the 2015 adults, for example, OK was the best model for both statistics but did not feature among the best for either of the two other years. Geostatistical methods, particularly OK, perform relatively poorly with clustered samples and high CoV (Webster and Oliver 2007). The clustered sampling of the 2007 data and the relatively high CoV for both the 2007 and 2015 datasets could explain the model performances in these years.

The inclusion of globally inputted environmental layers as covariates served to improve model performance but also likely added bias to the analysis. The inclusion of depth, habitat, and geographic location layers in the spatial analysis was justified based on the biology and ecology of the species, but these layers were also included because they were easy to obtain. The resolution of these layers, some of which were created for different purposes, did not necessarily reflect the scale at which ecological processes that affect queen conch distribution occur. As a result, undue bias may have been introduced, possibly affecting local autocorrelation and resulting in the false regression significance of some variables (Segurado et al. 2006). The effect of these biases is reflected in the output of the species distribution maps, on which some distributions appear to closely follow the trend of one or more of the covariates, so they poorly captured the local-scale spatial population structures. In addition, the significance of the covariates differed among the various models, as would be expected for models with different assumptions and specifications (Potts and Elith 2006). Therefore, model performance, and by extension the reliability of the estimated abundance index, may have been affected.

In the future, estimations could be further improved by obtaining geospatial layers of variables that reflect the resolution and scale of the biological processes of the species. Incorrect model calibration, unmodelled interactions and missing variables are also important sources of bias in spatial modelling (Segurado et al. 2006, Yu et al. 2013). Given the level of density-dependence and high exploitation of the queen conch (Stoner and Ray-Culp 2000), covariates explaining the impact of density-dependent mechanisms (such as the allee effect) and fishing on the species distribution may be important missing variables that could have improved the models. The systematic targeting of high-density patches of sedentary species is a key factor in modifying their spatial distribution directly by removing individuals and indirectly by affecting their habitat (Stoner et al. 2018). The inclusion of these two covariates would be prudent for future studies of this species and similar exploited sedentary species using regression-type spatial models.

Despite the improved access and availability of spatial interpolation techniques, there are fisheries such as that of the queen conch in which design methods of abundance indices estimation are widely used. Queen conch is a CITES (Convention on International Trade in Endangered Species of Wild Fauna and Flora) Appendix II listed species, a conservation status which requires careful management if the fishery is to be sustained (Prada et al. 2017). Improving abundance index estimations and understanding the spatial structuring of juvenile and adult aggregations is therefore an important part of this management process and can be used to develop integrated ecosystem management.

An ecosystem approach to fisheries (Garcia et al. 2003) implies explicit inclusion of the spatial aspect of species management, including small and larger-scale spatial trends caused by biological, ecological and environmental factors. The position and extent of the aggregations realized by the distribution map isolines are indicative of habitat preferences of juveniles and adult spawning stocks that are critical for management (Fig. 6). The changes in size and structure of these aggregations over time are indicative of the distributional response of conch biomass to environmental changes and the level of exploitation. High-density adult aggregations from our models can also serve as a proxy indicator of fisher behaviour since fishers are likely to target these areas. The spatial structures produced from our analysis are therefore a means for prioritizing areas for spatial planning measures such as zonation and targeted research, for example, to maintain viable spawning stock biomass above a reference point. This type of focused spatial management can help to understand critical habitat requirements and maintain population structure and genetic connectivity between local and regional populations (Kitson et al. 2018).