Supporting Spatial Management of Data-Poor, Small-Scale Fisheries With a Bayesian Approach

Rehren, Jennifer; Pennino, Maria Grazia; Coll, Marta; Jiddawi, Narriman; Muhando, Christopher

doi:10.3389/fmars.2021.621961

ORIGINAL RESEARCH article

Front. Mar. Sci., 01 July 2021

Sec. Marine Fisheries, Aquaculture and Living Resources

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.621961

This article is part of the Research Topic Using Ecological Models to Support and Shape Environmental Policy Decisions View all 28 articles

Supporting Spatial Management of Data-Poor, Small-Scale Fisheries With a Bayesian Approach

$\r\nJennifer Rehren*$ Jennifer Rehren^1*

Maria Grazia Pennino¹

Marta Coll²

Narriman Jiddawi³

Christopher Muhando⁴

¹Centro Oceanográfico de Vigo, Instituto Español de Oceanografía (IEO), Vigo, Spain
²Institute of Marine Science, Spanish National Research Council (ICM-CSIC), Barcelona, Spain
³Institute of Fisheries Research MALNF, Zanzibar, Tanzania
⁴Institute of Marine Sciences, University of Dar es Salaam, Zanzibar, Tanzania

Marine conservation areas are an important tool for the sustainable management of multispecies, small-scale fisheries. Effective spatial management requires a proper understanding of the spatial distribution of target species and the identification of its environmental drivers. Small-scale fisheries, however, often face scarcity and low-quality of data. In these situations, approaches for the prioritization of conservation areas need to deal with scattered, biased, and short-term information and ideally should quantify data- and model-specific uncertainties for a better understanding of the risks related to management interventions. We used a Bayesian hierarchical species distribution modeling approach on annual landing data of the heavily exploited, small-scale, and data-poor fishery of Chwaka Bay (Zanzibar) in the Western Indian Ocean to understand the distribution of the key target species and identify potential areas for conservation. Few commonalities were found in the set of important habitat and environmental drivers among species, but temperature, depth, and seagrass cover affected the spatial distribution of three of the six analyzed species. A comparison of our results with information from ecological studies suggests that our approach predicts the distribution of the analyzed species reasonably well. Furthermore, the two main common areas of high relative abundance identified in our study have been previously suggested by the local fisher as important areas for spatial conservation. By using short-term, catch per unit of effort data in a Bayesian hierarchical framework, we quantify the associated uncertainties while accounting for spatial dependencies. More importantly, the use of accessible and interpretable tools, such as the here created spatial maps, can frame a better understanding of spatio-temporal management for local fishers. Our approach, thus, supports the operability of spatial management in small-scale fisheries suffering from a general lack of long-term fisheries information and fisheries independent data.

Introduction

Small-scale fisheries employ over 90% of the world’s capture fishers (FAO., 2015, 2018) and are the major livelihood and protein suppliers in many coastal communities around the world (Chuenpagdee, 2011; Belhabib et al., 2015; Teh and Pauly, 2018; Loring et al., 2019; Salas et al., 2019). It is believed that well-managed small-scale fisheries can contribute to poverty alleviation and food security (Bene et al., 2007; Purcell and Pomeroy, 2015). However, assessing and managing these fisheries is challenging given the large number of species caught and the adaptive behavior of fishers in space, time, and fishing methods (Wiyono et al., 2006; Salas et al., 2007; Daw, 2008). The lack of alternative livelihoods and the strong resource dependency of many small-scale fishing communities impede common management measures such as total allowable catches or effort regulations (Pomeroy, 2012). Within the context of a global agenda to protect 10% of coastal and marine ecosystems through area management by 2020 (CBD, 2010), many tropical countries attempt to manage their coastal areas through different use-zones (Wells et al., 2007; De Santo, 2013).

Such an example is found in Zanzibar (Tanzania), where most of the coastline has been designated a conservation area ranging from general use zones to locally managed partially protected and privately managed no-take areas (McLean et al., 2012; Rocliffe et al., 2014). Zanzibar has achieved international targets by protecting 11% of its continental shelf, but a rapid appraisal by regional experts estimated that only 25% of the coral reef MPAs are effective (Rocliffe et al., 2014). Chwaka Bay on the east coast of Zanzibar is an important, year-round fishing area, which is part of Zanzibar’s large Mnemba Island Marine Conservation Area management plan (MIMCA) (McLean et al., 2012). But compliance with mesh-size and gear regulations is low (de la Torre-Castro and Lindström, 2010; Wallner-Hahn et al., 2016), making the bay a general use zone. A long history of intense exploitation (de la Torre-Castro and Rönnbäck, 2004; Rehren et al., 2018a), an increase in fishing effort (de la Torre-Castro and Lindström, 2010; Department of Fisheries Development., 2016), the use of illegal gears, and spatial use-conflicts (de la Torre-Castro and Lindström, 2010) have led to concerns for the sustainability of Chwaka Bay’s fisheries. In a participatory workshop in 2016, invited fishers advocated for implementing a no-take zone to combat the decrease in their catches and the reoccurring user conflicts (Rehren, 2017).

However, a prerequisite for the success of such no-take zones is to understand the spatial distribution of target species and identify its environmental drivers. While the people of Chwaka Bay strongly depend on fisheries resources for livelihoods and food security (Jiddawi, 2012), fisheries managers face scarcity and low-quality of data (Rehren et al., 2020). Because of the high-cost and spatial limitations of fisheries independent data collection, often the only source of information is landings data of individual fishers. This information is relatively easy to collect but comes with a strong sampling bias (Pennino et al., 2019). Spatio-temporal modeling approaches, therefore, need to account for all dependencies in the data, use information from different sources, and quantify associated uncertainties. The latter is particularly important to better understand the risks related to management interventions. Bayesian hierarchical species distribution models are well suited for this purpose because they allow for a more accurate estimation of uncertainty, given that observed data and model parameters can be considered as random variables (Banerjee et al., 2004).

We use a Bayesian hierarchical species distribution modeling approach on landing data from different fishing gears collected in 2014 to assess and predict the distribution range of key target species of Chwaka Bay. We identify common environmental drivers of distribution and areas of overlapping high relative abundance to prioritize potential conservation areas. The analyzed species represent key target resources found in fisheries catches throughout the Western Indian Ocean. Thus, our results serve as baseline information for future studies in the region.

Materials and Methods

Study Area

Chwaka Bay is a semi-enclosed bay-system located on the East Coast of Zanzibar (Tanzania) (Figure 1). The bay is relatively shallow, with depths up to 20 m in the outer borders and some parts of the bay falling dry during low tide. The sea surface temperature ranges from 25 to 31°C and salinity from 35‰ at the bay opening to 26‰ in the bay proper (Jiddawi and Lindström, 2012). Strong tidal currents, with a mean tidal range of 3.2 m (Nyandwi and Mwaipopo, 2000), cause high turbidity in the bay by stirring up sediments (Gullström et al., 2006). The north-eastern (November–March) and south-eastern (April–October) monsoons drive the bay’s climate, with the latter showing stronger winds, longer rain periods, and lower temperatures (Shaghude et al., 2012). The bay consists of a large mangrove forest on the southern shore, dense seagrass meadows throughout the bay, and a fringing reef at the bay opening. These habitats form a continuum through particulate organic matter exchange (Mohammed et al., 2001) and tidal, seasonal, foraging, and ontogenetic migration of fish (Gullström et al., 2012).

FIGURE 1

Figure 1. Chwaka Bay, Zanzibar (Tanzania). The bay comprises large mangrove stands in the south, a fringing reef at the bay opening, and coral patches inside the bay. Seagrass meadows are found throughout the bay with dense aggregations toward the central part.

The diversity of habitats and the protection from wave energy through the fringing reef give rise to a highly productive, year-round fishing area surrounded by several fishing villages (Figure 1). The local community highly depends on the fisheries’ resources for income and protein supply (Jiddawi and Lindström, 2012). The fishery targets multiple species ranging from invertebrates (e.g., sea cucumber, octopus), reef- and seagrass-associated fish (e.g., parrotfish and rabbitfish) to large pelagic species (e.g., mackerels and jacks). The main fishing gears are basket traps, dragnets, handlines, spears, and, to a minor extent, floatnets, longlines, fences, and gillnets (Rehren et al., 2018a). Dragnet fishers are mainly from Chwaka village located in the south of the bay, and their numbers have increased over the years (de la Torre-Castro and Lindström, 2010). The nets are weighted down with stones and dragged over the seafloor. Spatial use-conflicts arise from dragnets’ damage to sensitive habitats and basket traps from other fishers (Jiddawi and Ohman, 2002; Mangi and Roberts, 2006; de la Torre-Castro and Lindström, 2010). Following a prohibition of dragnets in 2001, the fishing grounds off Marumbi village were demarcated with buoys to ensure the protection of Marumbi fishers from dragnet fishing (de la Torre-Castro and Lindström, 2010). Despite this locally enforced zone, all gears are deployed throughout the entire bay. For over 20 years, fishers report decreases in their catch rates (de la Torre-Castro and Rönnbäck, 2004; Geere, 2014), which, together with the use of small mesh sizes and destructive gears, has led to a general concern of overfishing in the bay (de la Torre-Castro and Lindström, 2010; de la Torre-Castro et al., 2014).

Data Collection

Fisheries data, habitat, and depth information were collected by the first author during the north–east monsoon (January–June) and the south–east monsoon (September–December) season in 2014. Data collection was carried out on 18 days per month at the main landing sites (i.e., Chwaka village, Uroa village, and Marumbi village, Figure 1), covering a minimum of 30% of the fishing boats that went fishing on the day of sampling. The number of fishers sampled per gear and landing site was based on the gear and landing site’s relative proportion. The catch was classified to family, or if possible to species level (Bianchi, 1985; Anam and Mostarda, 2012), weighed to the nearest 1 g, and standardized to weight per fisher [weight per unit of effort (WPUE), kg fisher^–1]. The data collection was done directly at the beach during landing before the fishers sold their catch. Individuals of any size caught during fishing were landed and used at least for home consumption. The number of fishers, boat, gear, and fishing hours and the type of gear, boat, and propulsion were also collected. We assigned each sample to the corresponding lunar cycle (i.e., full moon, third quarter, new moon, and first quarter) and season (i.e., north–east monsoon and south–east monsoon). Information about the fishing location was collected as the name of the fishing ground. In April and December 2014, the main fishing grounds (71%) were mapped together with experienced fishers. The depth and seagrass and sand percentage cover were collected for 57% of the mapped fishing grounds and on additional non-fished, random locations in the bay. Depth was measured with a diving computer and corrected with the tide level records obtained from the Tanzania Ports Authority.¹ The substrate percentage cover was estimated within 2–6 quadrats at each location. Depth, seagrass, and sand were then interpolated within the spatial extent of the sampling locations using kriging techniques. Seagrass nor sand displayed any trend, and thus ordinary kriging was used with a spherical and exponential variogram model, respectively. Depth distribution showed a clear trend, and thus universal kriging with a cubic variogram model was applied to data detrended by a second-order polynomial trend surface analysis. The interpolation was done using the geoR package (Ribeiro et al., 2020). Shapefiles of coral reef presence–absence were obtained from the Institute of Marine Science, Zanzibar. The daily sea-surface temperature of 2014 was obtained from the GHRSST level 4 data set (0.01° × 0.01°) downloaded from the OPeNDAP data repository² using the XML package (Temple Lang, 2020) implemented in the R software. The temperature was transformed from Kelvin to Celsius, and the annual average was calculated. We used the habitat variables, depth, temperature, lunar cycle, and season as explanatory variables for the distribution of the WPUE of the most dominant species in the catches (Table 1). We chose the above explanatory variables because they have been identified as key predictors to determine spatial patterns of marine species (e.g., Beger and Possingham, 2008; Moore et al., 2009; Roos et al., 2015). The spatial distribution of the WPUE values was used as a proxy for the species’ relative abundance.

TABLE 1

Table 1. Key characteristics of the analyzed target species.

Statistical Analysis

All analyses and graphics were performed in R (R Core Team, 2020³). Prior to the analysis, the explanatory variables were standardized (i.e., difference from the mean divided by the corresponding standard deviation) (Gelman et al., 2014) using the decostand function in the vegan package (Oksanen et al., 2019) to better interpret both the direction (positive or negative) and magnitude (effect sizes) of the parameter estimates. We used the variance inflation factor (VIF, <3) (Zuur et al., 2009) and the Pearson correlation statistic to exclude covariates with high multi-collinearity. Only sand and seagrass were correlated for all species, and thus sand was removed from the analysis. Categorical variables were examined for an imbalance in the number of observations. Potential interactions between the response and a predictor covariate conditioned on other covariates were examined using coplots. Interactions between covariates were included in the model selection process if clear changes in the slope were observed, and the number of observations in each group was good enough to allow for such an analysis.

The relationship between the logarithm of the WPUE and predictors was modeled using a normal distribution (Figure 2). We included an independent and identically distributed random metier effect (Gómez-Rubio, 2020) that accounts for variations in WPUE due to differences in fishing methods and technologies (hereafter metier effect). Metiers were assigned to the different samples based on the associated fishing village,⁴ vessel, gear, and propulsion type. We further accounted for spatial autocorrelation by including a numeric vector with a mean of 0 and a Matern covariance function linking each observation to a spatial location. Thus, our model accounts for independent, region-specific, and metier-specific noise not explained by the available covariates. For the parameters involved in the fixed effects, vague Gaussian priors with a mean of 0 and a variance of 100 were used, while a gamma prior distribution on the precision τ with parameters 1 and 0.00005 was used for the metier effect. The random spatial effect only depends on two hyperparameters: the range and the variance of the spatial effect. Penalized complexity priors (Fuglstad et al., 2018) were used to describe prior knowledge on these hyperparameters. We set a prior range of 1 km with a probability of 0.05 for it to be lower and a prior variance of 1.7–2 (depending on the species) with a probability of 0.05 for it to be higher. We performed a sensitivity analysis of the choice of priors for the spatial effect by testing different priors and verifying that the posterior distributions were consistent and concentrated well within the support of the priors (Supplementary Figure 1).

FIGURE 2

Figure 2. Graphical representation of the model. The logarithm of the WPUE (γ_i) follows a normal distribution. For the fixed effects parameters, vague Gaussian priors with a mean of 0 and a variance of 100 were used. ζ_i is a Gaussian distributed random metier effect with a mean of 0 and a precision τ_ζ. By default, INLA assigns a gamma prior with parameters 1 and 0.00005 to the precision. The random spatial effect (ω_i) is approximated with a Matern covariance function (Q). The parameters κ and τ_ω determine the range and the total variance of the spatial effect. The penalized complexity priors of these parameters follow a normal distribution. This adaption of a Kruschke style diagram was generated using Bååth (2013) template for LibreOffice.

Bayesian inference was performed using the Integrated Nested Laplace Approximations (INLA) approach (Rue et al., 2009) with its corresponding package.⁵ INLA uses the so-called Stochastic Partial Differential Equation approach to approximate the Gaussian field with the Matern covariance function by a Gaussian Markov random field (Rue et al., 2009).

We selected the most parsimonious model, starting with all covariates (except those with VIF values > 3), based on the goodness-of-fit using the deviance information criterion (DIC) (Spiegelhalter et al., 2002) and Watanabe–Akaike information criterion (WAIC) (Watanabe, 2010; Supplementary Table 1). The model selection process was automated by using the Bdiclcpomodel_stack function available on GitHub.⁶ We included covariates in the final model if the probability for the regression parameters’ posterior distribution to be below or above zero was 80% or higher (depending on the relationship).⁷ The final model was evaluated with the log-conditional predictive ordinate (log-CPO) (Roos and Held, 2011), which is a “leave-one-out” cross-validation index to assess the predictive power of the model (Pennino et al., 2019). The CPO values were also used to identify outliers. We further evaluated the final model through residual plots (homogeneity of variance, outliers) (Supplementary Figure 2) and visualizing model predictions. Model assumptions were also analyzed by visualizing the predictive integral transform (Supplementary Figure 3), which measures the probability of a new value to be lower than the observed value (Held et al., 2010). INLA has built-in functions allowing for a linear interpolation of the spatial effect within each triangle into a finer regular grid. The resulting high-resolution map of the spatial effect can be seen as a proxy for the species’ relative abundance.⁸ The spatial effect maps were then stacked, and the posterior distribution of the mean, first, and third quantile was calculated to identify overlapping areas of high relative abundance.

Figures were created using the ggplot2 package, and maps were created with the marmap (Pante et al., 2020), mapdata (Richard A. Becker and Ray Brownrigg, 2018), mapproj (R by Ray Brownrigg McIlroy et al., 2018), cowplot (Wilke, 2019), ggspatial (Dunnington, 2018), rnaturalearth (South, 2017a), and rnaturalearthdata (South, 2017b) packages.⁹

Results

Drivers and Distribution of Target Species

Spatial dependencies and the random metier effect contributed strongly to the explained variance and hence improved model performance. The temporal covariates had a relatively strong effect on species distribution. While the third quarter and full moon were important covariates for Siganus sutor and Lethrinus lentjan (Figure 3), season was only important in the distribution of Leptoscarus vaigiensis. This relationship was positive, indicating higher relative abundances for L. vaigiensis during the south–east monsoon season (Figure 3). The strong positive relationship between the third quarter moon phase and the distribution of L. lentjan indicates that this is a relevant predictor of high WPUE values. Important environmental drivers were highly variable among species, but the magnitude of their effects was relatively similar (Figure 3). While depth was an important predictor for the distribution of L. vaigiensis, L. lentjan, and Lutjanus fulviflamma, temperature was important for Scarus ghobban, L. vaigiensis, and Lethrinus mahsena (Figure 3). S. sutor showed a much smoother spatial trend in relative abundance across the bay than the other species and had higher numbers of observations (twofold) (Figure 4). This species is highly dominant in the catches throughout the year and was particularly caught north of Uroa. The spatial pattern of relative abundance generally shows a clear south to north trend (Figure 4). None of the environmental variables were found to be important in the distribution of S. sutor (Figure 3), and only full moon was selected as a relevant predictor forming a negative relationship with abundance. The emperor species showed a strong positive relationship with depth (Figure 5), leading to higher relative abundances around the bay opening close to Uroa and Michamvi (Figure 4). Reef occurrence was the strongest predictor of L. mahsena (Figure 3), intensifying the spatial pattern of high relative abundances around the bay opening where the fringing reef is present. The distribution of L. fulviflamma and S. ghobban, contrastingly, showed a spatial pattern of higher relative abundances in the south of the bay (Figure 4). Although the distribution of L. fulviflamma was driven by greater depth, it is even stronger driven by high seagrass cover (Figure 5). It is present in dense seagrass meadows in the bay proper and around Chwaka village. The relative abundance of L. vaigiensis indicated a slightly negative relationship with seagrass and a strong positive relationship with depth (Figure 5), leading to a spatial pattern of high values around Uroa toward the bay opening and in the middle of the bay (Figure 4). The posterior distribution of the standard deviation of all spatial effects was relatively high, with the highest values toward the bay opening and the south of the bay, where the number of observations for the target species decreases (Figure 4).

FIGURE 3

Figure 3. Summary of the selected environmental drivers for each species and the value of the corresponding slope. The legend represents the probability of the slope (Importance) to be below (negative, red) or above (positive, blue) zero.

FIGURE 4

Figure 4. Posterior mean and standard deviation of the spatial random effect. Letters indicate the position of the villages: C, Chwaka; M, Marumbi; U, Uroa; Mi, Michamvi. Black lines depict the reef.

FIGURE 5

Figure 5. Functional response of the weight per unit of effort of each species to their main environmental drivers. Solid lines and shaded regions are the mean and 95% credibility intervals, respectively. NEMS stands for North–East Monsoon Season and SEMS for South–East Monsoon Season.

Identifying Areas of High Relative Abundance

Two main areas of high relative abundance were found by overlaying the mean posterior distribution of the spatial effect from all analyzed species: one area in the north of Uroa village and one area in front of Marumbi village (Figure 6). High relative abundances of the target species were also found close to the patch reefs inside the bay, which is a fishing area frequently visited by fishers even under unfavorable wind conditions as it is partly protected by the fringing reef. The fringing reefs and the deeper outer parts do not seem to create areas of particular high relative abundance for the analyzed target species. The area in front of Marumbi corresponds to the demarcated dragnet-free zone enforced by the Marumbi villagers. The area in the north of Uroa lies in the region suggested by the fishers as a no-take zone during the participatory workshop conducted in 2016. Both also remain areas of higher relative abundance in the posterior distribution of the upper and lower quantiles.

FIGURE 6

Figure 6. Posterior distribution of the mean, first, and third quantile of the combined spatial effect of the analyzed species to identify areas of overlapping high relative abundance. The dots represent a dragnet free zone (blue) enforced by the Marumbi fishers and the area proposed for the implementation of a no-take zone during the participatory workshop in 2016 (black). Letters indicate the position of the villages: C, Chwaka; M, Marumbi; U, Uroa; Mi, Michamvi. Black lines depict the reef.

Discussion

Main Drivers of Species Distribution

Understanding species distribution is a key aspect in setting successful spatio-temporal management plans (Franklin, 2009; Lawler et al., 2011; Guisan et al., 2013). The results from this study provide information on important distribution drivers of the WIO’s key target species and commonalities among them.

The environmental drivers found to be important for species distribution were highly dissimilar between the different species. This was very apparent among species of the same family: for the two emperor species, relevant selected predictors were opposite and for the two parrotfish, the only common environmental driver was temperature. Contrastingly, season was only an important driver of the distribution of L. vaigiensis. The weak influence of seasons on fish density is a general pattern observed in Chwaka Bay’s mangrove creeks (Mwandya et al., 2010) and the seagrass meadows close to Chwaka and Marumbi village (Lugendo et al., 2007). Only the heavy rain season from April to May has been shown by Lugendo et al. (2007) to induce significant changes in environmental factors and fish density inside the mangrove creeks of Chwaka Bay. This suggests relatively stable annual catches of the other target species, representing a significant part of trap and dragnet catch (Rehren et al., 2018a). Reduced temporal fluctuations in catches and the protection from wave energy through the fringing reef makes the bay a vital fishing ground that decreases the vulnerability of the fishing community. Attempts to reallocate fishing efforts to offshore areas, which has been part of past management actions (Gustavsson et al., 2014), need to compensate for a potential increase in income uncertainty.

The full moon lunar phase was the only relevant predictor for S. sutor. This species is of high importance to Chwaka Bay’s fishery since it dominates the main gears, and its annual yield strongly exceeds all other target species (Rehren et al., 2018b). S. sutor grazes over algae beds, and juveniles mainly use seagrass meadows as nurseries (Dorenbosch et al., 2005; Lugendo et al., 2005; Kimirei et al., 2011). Larger individuals are usually found around reefs and are associated with deeper depths (Kimirei et al., 2011). Accordingly, the spatial distribution of S. sutor shows a clear increasing trend in relative abundance toward the bay opening and, thus, deeper areas. S. sutor, unlike the more sedentary emperor species, is a relatively mobile species with a home range of about 900 m (Ebrahim et al., 2020a). This might explain its smoother distribution and the lack of patches in the spatial random effect compared to the other species.

Our analysis indicates that seagrass plays an important role in the distribution of the emperor species L. mahsena, the parrotfish L. vaigiensis, and the snapper L. fulviflamma. L. fulviflamma uses seagrass meadows and particularly mangrove swamps as nursery areas (Lugendo et al., 2005; Kimirei et al., 2011), explaining the higher relative abundances found in the south of the bay. During the workshop in 2016 with 30 participants, fishers reported that L. fulviflamma used to occur in higher quantities in the bay and that the species seemed to have moved toward deeper waters due to an increase in water temperatures (Rehren, 2017). While our model indicates a positive relationship of L. fulviflamma distribution with depth, temperature was not selected as a relevant predictor. Along this line, other studies conducted in Tanzania mainland found that depth best explained the size-frequency distribution of L. fulviflamma among habitats (Kimirei et al., 2015) with adult specimens found on deeper reefs (Kimirei et al., 2011). Temperature, however, was selected as the main driver for three of the other analyzed species, including S. ghobban, and overall species density is also negatively related to temperature in mangrove and mud/sand habitats of the bay (Lugendo et al., 2007). Dorenbosch et al. (2005) found high juvenile densities (>70%) and intermediate adult densities (30–70%) of S. ghobban in seagrass meadows, likely explaining its high relative abundance found in the central bay and in the south of the bay where seagrass meadows occur.

Lethrinus mahsena, also found to be driven by seagrass cover, is a generalist occurring in all habitats of the bay (Dorenbosch et al., 2005) and is particularly associated with coral patches and fringing reefs adjacent to seagrass beds (Gell and Whittington, 2002; Locham et al., 2010). This observation also matches our findings that the most important driver of its abundance was reef, followed by temperature and seagrass. Accordingly, areas of high relative abundance of L. mahsena were found around coral patches inside the bay, which are surrounded by large seagrass meadows and at the fringing reef in the north of Uroa. High relative abundances were also found in front of Marumbi village, a fishing ground dominated by dense seagrass beds.

Little information was available for the other emperor species, L. lentjan, which occurs in all habitats of the bay (Lugendo et al., 2005). The adult part of the population mostly occurs around the reef areas (Dorenbosch et al., 2005). Depth was the only environmental predictor selected in our model and probably explains the higher relative abundances in the north of Uroa and the area around Michamvi.

While L. vaigiensis mainly occurs in seagrass beds (Dorenbosch et al., 2005; Lugendo et al., 2005) and feeds on seagrass plants (Gullström et al., 2011), our results showed a slight negative relationship between seagrass cover and relative abundance. Although these results seem counterintuitive, Gullström et al. (2011) also found a negative relationship between shoot density and the variability of juvenile and adult density of L. vaigiensis. Fish assemblages in coral reef and seagrass habitats in Kenya likewise showed a decrease in overall density with increasing seagrass density (Chirico et al., 2017). The authors argued that this relationship possibly arises due to the reduced movement ability of fish in very dense and relatively short seagrass beds. Stronger environmental drivers of L. vaigiensis were depth and temperature in our models, which probably explains its high relative abundance around the north of Uroa and in the middle of the bay. Gullström et al. (2011) also found temperature to be a driver for the abundance of L. vaigiensis within different seagrass meadows in Chwaka Bay, but not depth. The authors, however, mainly studied seagrass meadows along the shore, which does not represent the full range of depth in the bay, possibly explaining the differences in our model results.

Potential Areas for Conservation

Spatio-temporal management is a key strategy to help mitigate conflicts among fishers and protect essential habitats and target species, without entirely depriving fishers of their economic basis (Rassweiler et al., 2012; Kerwath et al., 2013; Sale et al., 2014; Di Franco et al., 2016; Sala et al., 2021). In the WIO region, the implementation of conservation areas has been a prime management tool to reduce anthropogenic pressures (IUCN., 2004). In this study, we identified areas of high relative abundance of six key target species of the region to provide information for the prioritization of such conservation areas.

The identified overlapping areas of high relative abundance are in the north of Uroa, close to reef areas, and in front of Marumbi village, dominated by seagrass meadows. Furthermore, areas close to the patch reefs surrounded by seagrass meadows inside the bay also showed higher relative abundances. Both emperor species, the rabbitfish S. sutor, and the parrotfish L. vaigiensis would benefit from the closure of fishing in the selected areas. Except for Uroa, the identified areas do not occur on the proper fringing reef that runs along the bay opening. Although we use the analyzed WPUE values as a proxy for relative abundance, the absence of high relative abundances on the fringing reef is likely a mere reflection of the distribution of effort: the exposure and deeper depths at the fringing reef make it harder to fish with the main fishing methods. However, the WPUE distribution pattern indicates that proper reef areas with high coral cover are not necessarily areas with the highest fishing pressure in small-scale fisheries of the WIO region and that non-reef areas in Chwaka Bay may be as suitable for spatio-temporal management plans. These findings match the observation from de la Torre-Castro et al. (2014) that seagrass meadows, and not reefs, are the fishing grounds with the highest community benefits for the fishers of Chwaka village. In the WIO region, spatial management plans are often implemented to protect a specific habitat (Turpie et al., 2000; Wells et al., 2007; Rocliffe et al., 2014), which has led to the disproportionate representation of coral reefs in marine conservation areas (Wells et al., 2007; de la Torre-Castro et al., 2014; Chirico et al., 2017). The need to include seagrass meadows into fisheries management efforts has not only been highlighted for the bay (de la Torre-Castro et al., 2014) but globally (Unsworth et al., 2019).

Conservation areas are often selected without incorporating fisher’s behavior in the implementation of spatio-temporal management plans, which has led to weak compliance (McClanahan et al., 2006; Rosendo et al., 2011) and reduced benefits for fishing communities (Benjaminsen and Bryceson, 2012). For instance, Marine parks in Kenya have been established with little consultation of fishing communities, and in Tanzania, examples of opposing the enforcement of existing conservation areas exist (Wells et al., 2007). This is surprising as fishers have shown to possess strong ecological knowledge about their target stocks (Silvano and Valbo-Jørgensen, 2008). Lopes et al. (2018) have shown that fisher’s knowledge can even be reliable enough for predicting species occurrence. These observations are also reflected in our analysis: the two areas of high relative abundance correspond to the areas that: (1) have been prioritized by fishers for the dragnet free zone in front of Marumbi village in 2001 (de la Torre-Castro and Lindström, 2010); and (2) have been proposed as a potential no-take zone in the workshop of 2016 (Rehren, 2017). This study has been conducted to support local spatio-temporal management actions with quantitative information. In a series of upcoming participatory workshops, the relative abundance maps with their associated uncertainties will be used to effectively visualize and convey our key findings to the local stakeholders. With these workshops, we aim to combine short-term fisheries dependent data and fishers’ knowledge to synthesize the most relevant information about the spatial dynamics of Chwaka Bay’s fisheries and target resources and to prioritize spatial management actions.

Potential and Limitations of the Approach

In many small-scale fisheries, dependence on resources for livelihood and protein supply is high, making their sustainable management particularly important (Belhabib et al., 2015; Teh and Pauly, 2018; Loring et al., 2019; Salas et al., 2019). Appropriate management plans are, however, impeded by the notorious lack of data (Salas et al., 2007; Salayo et al., 2008; Samoilys et al., 2015). Fisheries independent surveys are often cost-intensive and visual census data collected around Zanzibar are spatially limited and only represent a temporal snapshot (Rehren et al., 2020). Fisheries catch information, on the other hand, is collected throughout the WIO region at a subset of landing sites (UNEP-Nairobi Convention and WIOMSA, 2015) and thus becomes the most cost-effective and accessible source of information when evaluating the spatio-temporal dynamics of target resources. Suppose that fishers’ catches represent thousands of sampling observations (García-Quijano, 2007) and those fishers use multiple gears catching a multitude of species. In that case, it can be argued that fisheries catches as a whole might better reflect species relative abundance than spatially and temporally limited visual census data. Bayesian hierarchical modeling approaches can better handle problems associated with this data, such as spatial dependencies and the fisher effect, through their direct incorporation in the model formulation (Banerjee et al., 2004; Pennino et al., 2019).

Our analysis shows that a large part of the variance was explained by the random effect terms highlighting the importance of spatial dependencies and effects stemming from the use of different gears, boats, and propulsions. The latter effect is very high, suggesting that Chwaka Bay’s fisheries might be better managed based on fishing units. This requires flexible and adaptive management approaches tailored around the dynamic behavior of fishers in space, time, and fishing methods. A clear benefit of Bayesian models in data-poor situations is the provision of uncertainty associated with the data and the parameter estimates (Banerjee et al., 2004; Fonseca et al., 2017). This is particularly important when the stakes are high, as is the case in small-scale fishing communities. Our analysis shows that the uncertainties associated with our results are relatively high, particularly for the areas in the south and the north of the bay. A central issue associated with fisheries-dependent data is that fishers have prior knowledge of the probability of catching their target species at a given location leading to sampling bias. Furthermore, in our case study, greater depths and the presence of hard corals make fishing harder for dragnet fisher, lowering the number of observations at the fringing reef. Our approach does not account for such sampling bias in the data, which might have influenced the identification of the high relative abundance areas. Another limitation is the difficulty in obtaining a precise geo-localization of the catch in tropical, small-scale fisheries because of the large number of vessels, their dynamic behavior, and the lack of technical equipment. Usually, the spatial location of the catch is associated with a fishing ground name and mapped subsequently, or the fishing ground location is identified on a map by the fisher. These procedures reduce the spatial precision of the catch and can mislead inference. The observation that our model selects the same area for conservation as fishers did during the participatory workshop in 2016 (Rehren, 2017) increases the confidence in our model results. It must be noted that these action plans were formulated and discussed with a limited number of fishers. A comparison of our results with information from ecological studies about the habitat preference and ecology of four of the analyzed species also suggests that our approach predicts the distribution of the analyzed species reasonably well. For the remaining two species (i.e., S. ghobban, L. lentjan), not enough information was found to evaluate our models properly. For S. sutor, it is likely that we missed to include the distribution of macroalgae in the bay as a predictor variable because S. sutor grazes on epibenthic algae and feeds on macroalgae thalli (Ebrahim et al., 2020b). But it is also possible that the spatial change in environmental or habitat variables in the bay is not strong enough to significantly affect S. sutor’s distribution because the area is relatively small and S. sutor’s mobility is relatively high. Salinity and primary productivity are other predictors that have been identified to drive species distribution (Roos et al., 2015; Gonzáles-Andrés et al., 2016; Coll et al., 2019). Studies from the bay, however, suggest that salinity is not a driver of fish density (Lugendo et al., 2007) and that habitat variables generally are more important predictors for fish assemblages and abundance (Gullström et al., 2008; Mwandya et al., 2010).

Information on environmental drivers at a high-resolution scale is often lacking, making it difficult to model the distribution of resources in small fishing areas such as Chwaka Bay. We used data collected by the first author about depth and habitat variables and interpolated them to get an estimate at all fishing grounds. Thereby, we did not consider the uncertainty associated with the covariate estimations, which can cause erroneous inference and a biased estimate of the covariate effect (Martínez-Minaya et al., 2018). Consequently, we compared models with and without the interpolated data: while some species had similar spatial random effect maps, for other species using the non-interpolated data resulted in peaks or throughs at missing locations (Supplementary Figure 4). In other words, including only a subset of the data would have resulted in misidentifying areas of high relative abundance. Ideally, information on environmental covariates should be available at all fishing grounds to avoid potential misalignment.

Areas prioritized by the fishing community or the approach used here may not be sufficient to achieve ecological objectives. The structural complexity of seagrass meadows in the bay, for instance, is an important factor that can drive fish abundance (Gullström et al., 2008) and habitats often function together with surrounding habitats determining fish composition through seascape structure (Berkström et al., 2012). Furthermore, the areas prioritized by the Chwaka Bay fishers are relatively small, while large marine protected areas are associated with higher success, particularly when protecting highly mobile species (Claudet et al., 2008; Vandeperre et al., 2011; Edgar et al., 2014; White et al., 2017). But it has also been shown that small community-based marine protected areas established in coral reef developing nations may nonetheless be highly successful (Ban et al., 2011; Chirico et al., 2017). In the long-term such conservation efforts need to be scaled up to regional or national levels to achieve the ecological principles of complementarity, representativeness, and connectivity (Ban et al., 2011).

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/Jrehren/Frontiers-2020-Rehren-Supporting-spatial-management-with-Bayesian-approach.git.

Author Contributions

MP, MC, and JR conceived and designed the research. NJ and JR designed the data collection. MP and JR performed the analysis. CM provided data for the analysis and the maps. MP, MC, NJ, CM, and JR wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This project was funded by the German Research Foundation (DFG) within a research fellowship program (RE 4358/1-1). The data used in the analysis was collected in 2014 within a program (SUTAS) funded by the Leibniz Gemeinschaft. MC acknowledges the ‘Severo Ochoa Centre of Excellence’ accreditation (CEX2019-000928-S) to the Institute of Marine Science (ICM-CSIC).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.621961/full#supplementary-material

Footnotes

^ https://www.ports.go.tz/index.php/en/customer-center/sea-tide-tables?limitstart=0 (accessed 05.02.2015).
^ https://opendap.jpl.nasa.gov/opendap/
^ www.r-project.org
^ Potential factors associated with higher variability in WPUE between fishers of different villages includes the travel distance to less/more productive fishing grounds and the possession of different levels of knowledge about fishing grounds and species.
^ http://www.r-inla.org/
^ https://github.com/MgraziaPennino/SDMs-with-INLA
^ The only exception was the covariate seagrass in the model of Leptoscarus vaigiensis, which was selected by the goodness-of-fit indicators but had only a probability of 70% for the posterior distribution to be below zero. We kept this covariate given its importance as habitat and food item for Leptoscarus vaigiensis.
^ The data and the script for the model construction can be accessed at https://github.com/Jrehren/Frontiers-2020-Rehren-Supporting-spatial-management-with-Bayesian-approach
^ Source code for the map can be accessed at https://github.com/MgraziaPennino/Create_map_study_area/blob/master/Map_study_area.R

References

Anam, R., and Mostarda, E. (2012). Field Identification Guide to the Living Marine Resources of Kenya. FAO Species Identification Guide for Fishery Purposes. Rome: FAO.