| Citation: | Lisha Guan, Xianshi Jin, Tao Yang, Xiujuan Shan. Ensemble habitat suitability modeling of stomatopods with Oratosquilla oratoria as an example[J]. Acta Oceanologica Sinica, 2023, 42(4): 93-102. doi: 10.1007/s13131-022-2051-z
|
Stomatopods are better known as mantis shrimp and mostly inhabit shallow marine waters in temperate and tropical regions (Reaka, 1980). They are generally predatory crustaceans who compete with carnivorous fishes for a broad diet of prey, while supply food for many aquatic consumers (e.g., fishes, turtles and seabirds) (Yang, 2001; Ng, 2013; DeVries et al., 2016). Consequently, stomatopods play important roles in structuring the benthic communities in various marine ecosystems globally (Reaka et al., 2008; Antony et al., 2010; Gianguzza et al., 2019; Wu et al., 2019). Moreover, most stomatopods depend on specific sediment characteristics for the building and irrigation of their burrows, and can significantly affect marine nutrient cycling and energy flow through bioturbation of sediments (Covich et al., 1999; Laverock et al., 2011).
Some stomatopod species are not only important in ecology but exploited as valuable fishery resources, such as Squilla mantis and Erugosquilla massavensis in the Mediterranean (Maynou et al., 2004; Abo-Hashesh et al., 2020) and Oratosquilla oratoria in the Northwest Pacific (Nakajima et al., 2010; Kim et al., 2017; Zhang et al., 2018). Substantial fishing pressures have largely depleted the stock abundance and damaged the benthic habitats of such stomatopods, while they remain as a dominant species among various fisheries resources in many coastal waters worldwide (Kodama et al., 2006; Antony et al., 2010; Tao et al., 2018; Ragheb et al., 2019). In this context, scientific conservation of exploited stomatopods is essential for the sustainability of coastal fisheries. This requires accurate identification of their suitable habitats as a prerequisite, as exploited stomatopods mostly dwell in burrows and can affect the food web dynamics of local ecosystems (Antony et al., 2010; Gianguzza et al., 2019).
Oratosquilla oratoria is a typical species of commercially exploited stomatopods with a wide distribution in coastal waters of the Northwest Pacific (Manning, 1971). It supports important local fisheries in China, Japan and Korea (Nakajima et al., 2010; Kim et al., 2017; Zhang et al., 2018). Its dwelling habit and oceanic barriers to larval transport can lead to evident genetic differences among very local stocks of O. oratoria (Lui et al., 2010; Yang and Li, 2018). Such local adaptation may indicate spatially varying habitat preferences, given heterogeneity in environmental conditions and their combined non-linear effects on O. oratoria distribution and habitat suitability (Li et al., 2018, 2020). Yet, few studies have published to promote accurate habitat identification of O. oratoria and other stomatopods as well (Zhang et al., 2020), obstructing scientific management and conservation of these valuable organisms.
The aim of this study was to provide a robust modeling framework for habitat suitability modeling of O. oratoria and similar stomatopods at wide coastal waters, so as to support their conservation and related fisheries management. The O. oratoria stock in the Bohai Sea was chosen as an example to develop the modeling framework, because of high quality of available survey and environmental data. The modeling framework contains two major parts: (1) selecting environmental predictors (especially sediment-related metrics) that better characterize O. oratoria distribution and building separate habitat suitability models (HSM) utilizing two selected modeling techniques; (2) comparing performance of the individual HSMs and integrating them to check whether the ensemble model outperforms either individual model, according to fishers’ knowledge and scientific survey data. The model outputs were expected to support place-based management and conservation of O. oratoria stock in the Bohai Sea.
Oratosquilla oratoria has high local fishing importance in the Bohai Sea, China, due to the depletion of fish and invertebrate stocks of higher values. The landings peaked at about 220000 t in 2006 and then gradually decreased to about 100000 t in recent years. Monthly landings showed an obvious seasonality, with the highest market values in spring (especially in May) and no allowable catch (except those caught by monolayer gillnets and jigs) during the summer fishing moratorium of China from 2009 to 2020. This seasonal suspension of fishing activities has been extended to gillnets and prolonged by one month since 2017, starting from May 1 instead of June 1 in the Bohai Sea. This resulted in a big overlap between the major fishing season of O. oratoria and the summer fishing moratorium. Local fishers had persistently requested for fishing O. oratoria at the beginning of the moratorium period. Finally, the Ministry of Agriculture and Rural Affairs, China, issued special licenses for three place-based O. oratoria fisheries in the Bohai Sea from May 1 to May 15 in 2021.
These place-based O. oratoria fisheries kicked off spatial management and conservation of stomatopod stocks in China’s coastal waters. The open areas for these fisheries were mainly determined by local fishermen and fisheries managers, due to limited information available on the suitable habitats of O. oratoria. This could put the O. oratoria stock at risk of overexploitation and even challenge biological conservation of the whole Bohai Sea, as O. oratoria is the most dominant nektonic species and can easily affect the food web dynamics of the local ecosystem (Wu et al., 2019).
This study covered almost the whole Bohai Sea (Fig. 1). Some nearshore waters have been excluded, as they are often occupied by stationary fishing devices or not deep enough (<5 m) for bottom trawling of survey vessels. This coastal water has an average depth at about 18 m, and is characterized by low salinity (<32.5), eutrophication and soft bottoms (mainly silt, clay and sand).
Oratosquilla oratoria density data were obtained from monthly bottom trawl surveys conducted by the Yellow Sea Fisheries Research Institute in the Bohai Sea in 2017. These surveys followed a fixed-station sampling design, with one station at each intersection of grid of 0.25°×0.25° (longitude×latitude). Some grids within the study area were not sampled for various reasons. For instance, lots of scattered stationary fishing devices prohibited trawling along the mid-coast and outside the Bohai Bay especially in March and April, and some large marine debris aggregated in the inner Liaodong Bay, which could cause serious damages on trawl nets. In addition, the same pair of bottom trawlers were utilized in these surveys, following consistent sampling protocols. This pair of bottom trawlers are general fishing vessels with an engine power of around 147 kW and identical trawl nets of commercial size (i.e., a cod-end mesh size of 2 cm, a headline height of 5–6 m and a distance of 22.6 m between wings).
Only the survey data in spring (April to June) were utilized in this study, which include a total of 121 successful tows. Most of these tows continued for 60 min at an average speed of 3 kn. Exceptions occurred at some stations, where tow duration was reduced to about 30 min to avoid collisions with stationary fishing devices or moving boats. For each tow, catches of nektonic species were separated by species and individuals in each sample were counted and weighed. Catch per tow by species, in number and biomass, were recorded on the spot, along with the depth, geographic coordinates, date and time at the beginning and end of each tow. We calculated O. oratoria density for each tow by dividing the catch weight (g) of this species by tow duration (h).
The environmental predictors include depth, bottom temperature (BT), bottom salinity (BS), and six sediment content or grain-size variables. Depth, BT and BS were measured simultaneously with the bottom trawl surveys. For each tow, depths were measured with a sonar detector at the beginning and end of the tow, and then were averaged to represent the depth throughout the tow duration. Surface-to-bottom temperature and salinity profiles were obtained with a SEABIRD CTD in a few minutes after the bottom trawlers moving away from the end of each tow. BT was estimated as the average value of temperatures for the bottom layer of two-meters’ height, and BS was calculated likewise. Moreover, the ordinary Kriging method using the Gaussian semi-variogram function was applied to create bathymetry, monthly salinity and temperature fields for predicting and mapping habitat suitability index (HSI) of the study area from April to June, 2017. These fields were generally gridded by 0.05°×0.05° in latitude and longitude.
Two types of sediment metrics were considered in the study, which include three variables representing the percentages of sand, silt and clay (Psand, Psilt and Pclay for short) and three statistical parameters representing grain size distribution (namely mean grain size, skewness and kurtosis, i.e., Mz, Sk, and Kt) in the surface substrates (0–3 cm). One variable Psilt was ignored in this study because of its collinearity with other two variables Psand and Pclay.
Spatial distribution of the intervals of the five remaining sediment metrics were derived from recent literature on grain-size distribution of surface sediments in the Bohai Sea (Yuan et al., 2020). In this data source, Sk and Kt were calculated based on Mz and Mz was reported in Krumbein (1934) in phi (Φ) scale, where sand is Φ=−1 to Φ=4 (2 000 μm to 62.5 μm), silt is Φ=4 to Φ=8 (62.5 μm to 3.906 μm), and clay is Φ>8 (<3.906 μm) (Table 1). Consequently, larger Mz indicates smaller grain size (μm); positive and larger Sk indicates a longer tail in the right, representing smaller dominant grain size in Φ (i.e., larger dominant grain size (μm)) than mean grain size; and large Kt represents highly aggregated grain sizes around Mz.
| Grain size/mm | Grain size/μm | Wentworth class | Phi (Φ) range |
| 1 to 2 | 1000 to 2000 | very coarse sand | Φ=–1 to Φ=0 |
| 1/2 to 1 | 500 to 1000 | coarse sand | Φ=0 to Φ=1 |
| 1/4 to 1/2 | 250 to 500 | medium sand | Φ=1 to Φ=2 |
| 1/8 to 1/4 | 125 to 250 | fine sand | Φ=2 to Φ=3 |
| 1/16 to 1/8 | 62.5 to 125 | very fine sand | Φ=3 to Φ=4 |
| 1/32 to 1/16 | 31.25 to 62.5 | coarse silt | Φ=4 to Φ=5 |
| 1/64 to 1/32 | 15.625 to 31.25 | medium silt | Φ=5 to Φ=6 |
| 1/128 to 1/64 | 7.813 to 15.625 | fine silt | Φ=6 to Φ=7 |
| 1/256 to 1/128 | 3.906 to 7.813 | very fine silt | Φ=7 to Φ=8 |
| < 1/256 | < 3.906 | clay | Φ>8 |
DownLoad:
CSV
We replaced the five remaining sediment metrics with five integer variables of sediment (i.e., Ord_Psand, Ord_Pclay, Ord_Mz, Ord_Sk and Ord_Kt), each of which represented the orders of the intervals in ascending order for representative metrics (Table 2). These variables were assumed as static over short term (e.g., three to five years). We thereby derived five integer variables of sediment for sampling tows based on the midpoints of their beginning and end geographic coordinates.
| Ord_Psand | Sand content | Ord_Pclay | Clay content | Ord_Mz | Mean size (Φ) | Ord_Sk | Skewness | Ord_Kt | Kurtosis |
| 1 | 0 to 10% | 1 | 0 to 5% | 1 | Φ<3 | 1 | –2.5 to –1.5 | 1 | 0.5 to 1.5 |
| 2 | 10% to 20% | 2 | 5% to 10% | 2 | Φ=3 to Φ=4 | 2 | –1.5 to –0.33 | 2 | 1.5 to 2.5 |
| 3 | 20% to 30% | 3 | 10% to 15% | 3 | Φ=4 to Φ=5 | 3 | –0.33 to 0.33 | 3 | 2.5 to 2.75 |
| 4 | 30% to 40% | 4 | 15% to 20% | 4 | Φ=5 to Φ=6 | 4 | 0.33 to 1.5 | 4 | 2.75 to 3 |
| 5 | 40% to 50% | 5 | 20% to 25% | 5 | Φ=6 to Φ=7 | 5 | 1.5 to 2.2 | 5 | 3 to 3.25 |
| 6 | 50% to 60% | 6 | 25% to 30% | 6 | Φ>7 | 6 | 2.2 to 3.5 | 6 | 3.25 to 3.5 |
| 7 | 60% to 70% | 7 | 30% to 35% | – | – | 7 | > 3.5 | 7 | 3.5 to 4.5 |
| 8 | 70% to 80% | 8 | 35% to 40% | – | – | – | – | 8 | 4.5 to 5.5 |
| − | − | − | − | − | − | − | − | 9 | >5.5 |
DownLoad:
CSV
Eight environmental predictors were considered as preliminary predictor variables for habitat modeling, which included depth, BT, BS, Ord_Psand, Ord_Pclay, Ord_Mz, Ord_Sk and Ord_Kt. Two modeling techniques were chosen and applied for predictor selection and habitat suitability modeling of the targeted O. oratoria stock. One is generalized addition model (GAM), which can fit either presence-absence or density data and has shown good predictive performance on interpolation accuracy and transferability (extrapolative accuracy) in modeling species distribution or ecological niche (Heikkinen et al., 2012; Qiao et al., 2019). The other one is geographical weighted regression (GWR), which can estimate spatial non-stationary effects of environmental predictors on the distribution of dwelling crustaceans including O. oratoria (Li et al., 2018, 2020). GWR showed better interpolation accuracy than GAM in predicting the distribution of O. oratoria within a small coastal water in China (Li et al., 2020).
These two modeling techniques were firstly applied to determine which type of sediment metrics better characterize O. oratoria distribution and habitat and which predictors would be included for habitat suitability modeling based on their relative effects on O. oratoria distribution. Only the remaining predictors in the most parsimonious (or optimal) models were finally kept. The density of O. oratoria (D) was modeled as a function of predictor variables with a Tweedie error distribution and a log link function in GAM:
| $$ {g}\left({D}_{i}\right)={{{\text{λ}} }}_{0}+\sum _{k=1}^{m}{{s}}\left({x}_{ik}\right)+{\varepsilon }_{i}, $$ | (1) |
where xik is the value of the kth independent variable for tow i; m is the number of independent variables; λ0 is the intercept parameter; g() represents the link function; s() represents a spline smooth function for each numerical variable; and εi is the random error for tow i.
The densities of O. oratoria were scaled first and then modeled as a function of environmental predictors in GWR,
| $$ {\rm{lg}}\, D=\mathrm{lg}\left(D+1\right), $$ | (2) |
| $$ {{\rm{lg}}\, D}_{j}={\beta }_{j0}+\sum _{k=1}^{m}{\beta }_{jk}{x}_{jk}+{\varepsilon }_{j}, $$ | (3) |
where xjk is the value of the kth independent variable at location j; m is the number of independent variables as above; βj0 is the intercept parameter at location j; βjk is the local regression coefficient for the kth independent variable at location j; and εj is the random error at location j. A Gaussian kernel function was utilized in this model, with an optimal adaptive bandwidth (i.e., a fixed quantity that reflects local sample size) determined by an automatic bandwidth selection approach (Gollini et al., 2015).
A pseudo stepwise procedure was applied to find a subset of predictor variables for the most optimal models of GAM and GWR, respectively. The procedure proceeded in a forward direction following five steps (Lu et al., 2014): (1) start with bivariate GAM or GWR by sequentially regressing a single predictor variable against the dependent variable; (2) find the best performing model which produces the lowest Akaike Information Criterion (AIC) corrected for small sample size (AICc), and permanently include the corresponding predictor variable in subsequent models; (3) sequentially introduce a variable from the remaining group of predictor variables to build new models with the permanently included predictor variable, and identify the next permanently included variable from the best performing model with the lowest AICc; (4) repeat step (3) until all predictor variables are included; (5) find the optimal model with the lowest AICc.
We assumed a positive linear relationship between HSI and lg D for habitat suitability modeling, as the densities of O. oratoria had a highly right skewed distribution due to highly heterogeneous distribution of this specie. This could reduce the impact of large density values on underestimating HSI for most of the study area (Xue et al., 2017). A bivariate GAM was built to estimate the relationship between lg D and each environmental predictor included in the optimal GAM with the spline smooth regression. The estimated smooth regression was then used to establish suitability index (SI) curves for each predictor following the equation below:
| $$ {\rm{SI}}=\frac{\widehat{Y}-{\widehat{Y}}_{\min}}{{\widehat{Y}}_{\max}-{\widehat{Y}}_{\min}}, $$ | (4) |
where
| $$ {\rm{HSI}}=\frac{\displaystyle\sum _{k=1}^{n}{\mathrm{SI}}_{k}}{n},\quad k=1,2,\cdots ,n, $$ | (5) |
where n is the number of environmental predictors included in the optimal GAM. We called this model as GAM-based HSM hereinafter.
Additionally, we applied GWR to establish a set of SI curves for environmental predictors included in the optimal GWR model following Eq. (4), where
The GAM- and GWR-based HSMs were then applied separately to predict HSIs for the study area with the estimated or derived data of environmental predictors by grid within the domain. The spatial distributions of HSIs were mapped by month to observe the overall spatial trend in suitable habitat distribution for O. oratoria in spring. Interpolation accuracy of these HSMs were compared by regressing their respective predicted HSIs against O. oratoria densities at observed stations by month. Transferability of these models were evaluated by comparing HSIs to fishers’ knowledge at scarcely sampled areas within the study domain. Finally, HSIs from the GAM- and GWR-based HSMs were combined with the arithmetic mean method (Georgian et al., 2019) to show the effect of an ensemble model of these two models on suitable habitat identification.
The optimal GAM included five environmental predictors (i.e., BT, Ord_Sk, depth, BS and Ord_Mz) with the lowest AIC at 2203.717 (Fig. 2). All these predictors were statistically significant, with depth having p-value of its smooth terms equal 0.025 and four other predictors having p-values of their smooth terms smaller than 0.001. This optimal GAM explained 54.2% of the deviance in the O. oratoria density data. The optimal GWR model included four environmental predictors (i.e., BT, Ord_kt, depth and Ord_Mz), which had the lowest AICc at 538.656 (Fig. 2). The inclusion of additional predictors resulted in obvious increases in AICc (Fig. 2). This optimal GWR model explained 50.4% of the deviance in the scaled O. oratoria density data. Overall, the two sediment content predictors (Ord_Psand and Ord_Pclay) were both excluded in the optimal GAM and GWR model. In other words, grain-size metrics of sediment turned to be more important predictors for O. oratoria density distribution than the sediment content metrics.
The global effects (i.e., consistent effects over the study area) of included predictors on O. oratoria density displayed as the response curves from the bivariate GAM for each environmental predictor included in the optimal GAM (Fig. 3). The response curve of BT approximated an asymptotic curve with obviously high densities at temperatures between 13℃ and 22℃. The Ord_Sk response curve indicated greater O. oratoria densities with relatively larger (but not too large) Sk of the surface sediments. In other words, higher O. oratoria density appeared at sediments with relatively larger dominant grain size than mean grain size (μm). The response curve of depth suggested generally high O. oratoria densities within 25 m in depth. Additionally, the response curve of BS suggested that O. oratoria density increased fast with increasing BS before reaching 31 and then kept at high levels between 31 and 32.5. As for predictor Ord_Mz, its response curve indicated that O. oratoria density appeared the highest when Mz fell in the 4th interval (i.e., Φ=5 to Φ=6) with secondary higher densities when Mz fell in the 3rd and 5th intervals (i.e., Φ=4 to Φ=5 and Φ=6 to Φ=7; Yuan et al., 2020). This suggested that greater O. oratoria densities were associated at fine sediments (e.g., silt and mud with Mz between (Φ=4) and (Φ=7)) instead of coarse-grained sediments (e.g., sand with Mz<(Φ=4)) or much finer-grained sediments (e.g., clay with Mz>(Φ=8)).
Maps of the local coefficient estimates for environmental predictors from bivariate GWRs reflected their non-stationary effects on the density distribution of O. oratoria (Fig. 4). Positive relationships existed between O. oratoria density and BT, with the coefficients varying from 0.50 to 1.48 (Fig. 4a). The coefficients of BT were the highest in the northeastern Laizhou Bay and generally higher in the Liaodong Bay and the mouth of the Bohai Sea. The second included environmental predictor was Ord_Kt, which displayed strong positive effects on O. oratoria density in two coastal bays (i.e., the Laizhou Bay and the Liaodong Bay) and strong negative effects in the middle of the Bohai Sea (Fig. 4b). Depth showed strong positive effects on O. oratoria density in the Laizhou Bay and the Bohai Bay, while displayed generally negative effects in northeastern Bohai Sea where waters are mostly deeper than 25 m (Figs 1 and 4c). The final included predictor Ord_Mz displayed strong positively effects on O. oratoria density in the Liaodong Bay. Oppositely, it displayed strong negative effects in the inner Laizhou Bay and the middle Bohai Sea (Fig. 4d).
The GAM-based HSI predictions were the highest along the mid-coast of the Bohai Sea and secondary higher around the mouth of the Laizhou Bay (Figs 5a, c, e). Designed stations falling in these areas were seldom sampled due to widely dispersed stationary fishing devices (e.g., gillnets and trap nets). These areas matched well with the open areas for licensed mantis shrimp fisheries of adjacent provinces in May during the summer fishing moratorium of 2021 (Fig. 5c). In addition, habitat suitability tended to increase from April to June with the median HSIs by month rising from 0.53 to 0.59 to 0.63.
The GWR-based HSM displayed similar temporal pattern over months, with 0.56, 0.72 and 0.78 as medians in April, May and June. Differently, the GWR-based HSM identified the highest habitat suitability in the northeastern Laizhou Bay and on the east of the Bohai Bay (Figs 5d, f). Overall, the observed density distributions of O. oratoria appeared to match better with the GWR-based HSIs (Fig. 5). The linear relationships were statistically significant between log-transformed O. oratoria density and both model-derived HSIs. The relationship for GAM-based HSIs had little explanatory power (R2 = 0.064) on O. oratoria density distribution, whereas the linear relationship for GWR-based HSIs explained 25.2% of the variation in the observed O. oratoria densities (Fig. 6). On balance, integrating these two HSMs seemed to improve the predictive performance overall, successfully identifying fisher-recognized and survey-indicated suitable habitats for O. oratoria in either sparsely sampled or well investigated areas (Fig. 7).
The distribution of O. oratoria has been often associated with BT, BS, depth and metrics of sediment in the Yellow Sea and the Bohai Sea (Li et al., 2020; Zhang et al., 2020). This study identified BT, BS, depth and grain-size metrics of sediment as informatic predictors for habitat suitability of O. oratoria, while excluding the two metrics of sediment content (i.e., Ord_Psand and Ord_Pclay). This suggested that grain-size metrics of sediment might be better predictors than sediment content metrics for modeling the distribution or habitat suitability of sedentary species, especially those of burrowing habit (e.g., lobster, crab and mantis shrimp). For one reason, grain-size metrics reflect the hardness and texture of substrates, which affect the burrow development of cave dwellers (Chen et al., 2017). For another reason, sediment grain-size metrics can influence the thermoregulation of cave dwellers. Just as very loose or dense soil is not good for water storage (Zhao et al., 2019), grain-size metrics to a great extent determine soil water storage capacity. Soil heat capacity is commonly estimated as a function of the ratios of inorganic solid, organic solid, liquid and gas (Chen et al., 2019). Liquid has higher specific heat capacity than organic solid, inorganic solid and gas. Therefore, grain-size metrics directly link to the ratio of liquid contained and can affect sediment heat capacity, thereby influencing the thermoregulation of cave dwellers.
Different types of sediment metrics have been utilized to characterize the habitats of crustacean species of burrowing habit (e.g., lobsters, mantis shrimps and fiddler crabs). For example, one categorical variable of multiple sediment types (Tanaka and Chen, 2015), several numerical variables representing the contents of separate sediment contents (e.g., rock, silt or clay) (Li et al., 2020), and statistical metrics representing grain size distribution in the substrates (Li et al., 2018; Chen et al., 2019). Yet, related studies rarely reported any selection process among different types of sediment metrics. This study indicated that grain-size metrics outperformed numerical sediment content metrics on modeling habitat of O. oratoria. This suggested potential requirement to select optimal sediment metrics for improving habitat modeling of such cave-dwelling species.
The global effects of included predictors based on GAM are consistent with the local effects of corresponding predictors determined by GWR in most of the study area (Figs 3 and 4). First, increases in BT had generally positive effects on O. oratoria density at local scale. These positive effects appeared relatively weaker outside the Huanghe River Estuary, where BTs were obviously higher than other areas. This coincides the global effect of BT (i.e., O. oratoria density firstly increased and then kept relatively stable with increasing BT from 5℃ to 22℃). Likewise, the global and local effects of depth match well on O. oratoria density, suggesting that waters shallower than 25 m tended to be more suitable for this species than deeper waters. Zhang et al. (2020) reported a similar result in the Haizhou Bay (a large basin within the Yellow Sea) and adjacent waters, where the suitability of O. oratoria decreased with increasing depth.
In addition, the global effects of grain size metrics suggested that O. oratoria prefers sediments with mean grain sizes around the 4th level of Ord_Mz (i.e., silt and mud, Fig. 3; Yuan et al., 2020), and its sediment suitability would further increase with relatively larger dominant grain size than mean grain size (μm). These global effects of grain size metrics are consistent with their local effects. At local scale, higher density was linked to Ord_Mz at lower level (i.e., larger grain size (μm)) in the inner Laizhou Bay and the middle Bohai Sea, where Ord_Mz was mostly in higher levels (Levels 4 and 5) with finer grain size. Oppositely, higher density was associated with Ord_Mz at higher level (i.e., smaller grain size (μm)) in the Liaodong Bay with coarse-grained sediments and Ord_Mz in low levels (Levels 2−4) (Fig. 4d; Yuan et al., 2020). Overall, these local effects of Ord_Mz suggested favorable sediments of mean grain sizes around the 4th level of Ord_Mz (i.e., silt and mud).
Furthermore, larger skewness in grain size distribution (i.e., relatively larger dominant grain size than mean grain size) contributed to higher O. oratoria density at global scale. Differently, the kurtosis of sediment grain size distribution was identified as more important predictor than mean grain size for O. oratoria distribution at local scale. It requires further exploration to understand how kurtosis of sediment grain sizes could affect the distribution of O. oratoria in ecology. Furthermore, BS was included in the optimal GAM, but excluded in the optimal GWR model. Some recent studies (Li et al., 2020; Zhang et al., 2020) excluded BS as a significant predictor in modeling O. oratoria distribution or suitability in the Haizhou Bay and adjacent coastal waters, which was explained by this species’ ability to survive well in a wide range of salinity. Yet, the utilization of RNA-Sequencing techniques suggested that the most suitable salinity for O. oratoria cultured at 25℃ may be around 35 from the perspective of osmoregulation (Lou et al., 2019). In this context, the global effect of BS identified in this study did not contradict existing knowledge on O. oratoria ecology, given that BT ranged from 5−22℃ during the survey. The optimal GWR model excluding BS as an important predictor might be explained by low salinity variation at local scale.
The GAM-based HSM displayed better transferability with higher predictive accuracy in areas that were sparsely sampled, when compared to the GWR-based HSM. Suitable habitats predicted by the former model greatly overlapped with the open areas for licensed mantis shrimp fisheries during the summer fishing moratorium of 2021 (Fig. 5). These areas should support high abundance of mantis shrimp, as they were mainly determined by local fishermen and fisheries managers who closely track hotspots of their targeted species (Wilson, 2017). This inference was also supported by the fact that sampling in these areas were mostly prohibited by widely dispersed stationary fishing devices (e.g., gillnets and trap nets) in spring.
In addition, the GWR-based HSM outperformed the GAM-based HSM in interpolation accuracy, effectively quantifying habitat suitability in well sampled areas. High densities of O. oratoria were mostly observed in locations with high GWR-based HSIs. These locations mostly fall in the maritime jurisdictional zone of Shandong Province (Ding et al., 2020). Unfortunately, the application was not approved for licensed mantis shrimp fisheries in coastal areas along Shandong Province in 2021.
Improvement in the interpolation accuracy of the GWR-based HSM may be explained by well quantified spatial non-stationary habitat preferences due to local adaptation of O. oratoria. There are about 40 rivers entering into the Bohai Sea, forming various estuaries distributed well along the coast. Freshwater discharges sourced from different rivers usually vary in the amount of flux, contents of nutrients and sediment load. Weak hydrological connectivity between different estuaries and adjacent waters has led to highly heterogeneous environmental conditions in this ecosystem (Luo et al., 2021). This could have provided oceanic barriers to larval dispersal and promote local adaptation of O. oratoria, given its sedentary habit in burrows after the settlement of post-larvae (Ng, 2013). Moreover, local adaptation of O. oratoria was supported by the genetic diversity and differentiation in individuals collected in a local coast within the Bohai Sea (Liu et al., 2009), as well as the morphological differences between local stocks in the Laizhou Bay (a coastal bay within the Bohai Sea) and adjacent waters (Wang et al., 2020).
Furthermore, ensemble modeling technique enables a more robust characterization of ecological niche and habitat suitability by combining multiple types of HSMs to avoid dependence on a single model type (Robert et al., 2016). This explains why integrating the two individual HSMs successfully identified fisher-recognized and survey-indicated suitable habitats for O. oratoria in either sparsely sampled or well investigated areas. Recently, HSMs are getting wider applications into habitat evaluation and conservation of aquatic wildlife, as well as the forecasts of species invasions or distributional shifts in marine ecosystems. This has promoted increasing desire for higher accuracy in spatial predictions of habitat suitability, especially when occurrence or abundance data were sparse or unequally distributed over the study area of targeted species (Cerasoli et al., 2021; Qiao et al., 2019). In this context, ensemble modeling technique has been well utilized to identify suitable habitats for vulnerable corals in the South Pacific Ocean so as to support the design of new spatial closures for fisheries management in the region (Georgian et al., 2019).
In conclusion, this study provided an ensemble modeling framework for habitat suitability modeling of O. oratoria and similar stomatopods at wide coastal waters, so as to inform their conservation and related coastal fisheries management. This case study suggested that grain-size metrics of sediment may outperform sediment content metrics in modeling the habitat of cave-dwelling stomatopods. Moreover, the pros and cons of GAM- and GWR-based HSMs were illustrated and integrating them on balance could yield better performance in spatial predictions of suitable habitats at least for O. oratoria.
|
Abo-Hashesh A T, Madkour F F, Sallam W S, et al. 2020. Molecular identification and phylogenetic relationship of Erugosquilla massavensis (Kossmann, 1880) from the Mediterranean Sea, Egypt. Egyptian Journal of Aquatic Biology and Fisheries, 24(6): 533–545. doi: 10.21608/ejabf.2020.117576
|
|
Antony P J, Dhanya S, Lyla P S, et al. 2010. Ecological role of stomatopods (mantis shrimps) and potential impacts of trawling in a marine ecosystem of the southeast coast of India. Ecological Modelling, 221(21): 2604–2614. doi: 10.1016/j.ecolmodel.2010.07.017
|
|
Cerasoli F, Besnard A, Marchand M A, et al. 2021. Determinants of habitat suitability models transferability across geographically disjunct populations: Insights from Vipera ursinii ursinii. Ecology and Evolution, 11(9): 3991–4011. doi: 10.1002/ece3.7294
|
|
Chen Tung-Yun, Hwang Gwo-Wen, Mayfield A B, et al. 2017. The relationship between intertidal soil composition and fiddler crab burrow depth. Ecological Engineering, 100: 256–260. doi: 10.1016/j.ecoleng.2016.12.011
|
|
Chen Tung-Yun, Hwang Gwo-Wen, Mayfield A B, et al. 2019. The development of habitat suitability models for fiddler crabs residing in subtropical tidal flats. Ocean & Coastal Management, 182: 104931
|
|
Covich A P, Palmer M A, Crowl T A. 1999. The role of benthic invertebrate species in freshwater ecosystems: Zoobenthic species influence energy flows and nutrient cycling. Bioscience, 49(2): 119–127. doi: 10.2307/1313537
|
|
DeVries M S, Stock B C, Christy J H, et al. 2016. Specialized morphology corresponds to a generalist diet: linking form and function in smashing mantis shrimp crustaceans. Oecologia, 182(2): 429–442. doi: 10.1007/s00442-016-3667-5
|
|
Ding Qi, Shan Xiujuan, Jin Xianshi, et al. 2020. Research on utilization conflicts of fishery resources and catch allocation methods in the Bohai Sea, China. Fisheries Research, 225: 105477. doi: 10.1016/j.fishres.2019.105477
|
|
Franklin J. 2009. Mapping Species Distributions: Spatial Inference and Prediction. Cambridge: Cambridge University Press, 181–187
|
|
Georgian S E, Anderson O F, Rowden A A. 2019. Ensemble habitat suitability modeling of vulnerable marine ecosystem indicator taxa to inform deep-sea fisheries management in the South Pacific Ocean. Fisheries Research, 211: 256–274. doi: 10.1016/j.fishres.2018.11.020
|
|
Gianguzza P, Insacco G, Zava B, et al. 2019. Much can change in a year: The massawan mantis shrimp, Erugosquilla massavensis (Kossmann, 1880) in sicily, Italy. BioInvasions Records, 8(1): 108–112. doi: 10.3391/bir.2019.8.1.11
|
|
Gollini I, Lu Binbin, Charlton M, et al. 2015. GWmodel: An R package for exploring spatial heterogeneity using geographically weighted models. Journal of Statistical Software, 63(17): 1–50
|
|
Heikkinen R K, Marmion M, Luoto M. 2012. Does the interpolation accuracy of species distribution models come at the expense of transferability?. Ecography, 35(3): 276–288
|
|
Kim S, Kim H, Bae H, et al. 2017. Growth and reproduction of the Japanese mantis shrimp, Oratosquilla oratoria (De Haan 1844) in the coastal area of Tongyeong, Korea. Ocean Science Journal, 52(2): 257–265. doi: 10.1007/s12601-017-0027-2
|
|
Kodama K, Shimizu T, Yamakawa T, et al. 2006. Changes in reproductive patterns in relation to decline in stock abundance of the Japanese mantis shrimp Oratosquilla oratoria in Tokyo Bay. Fisheries Science, 72(3): 568–577. doi: 10.1111/j.1444-2906.2006.01185.x
|
|
Krumbein W C. 1934. Size Frequency distributions of sediments. Journal of Sedimentary Research, 4(2): 65–77
|
|
Laverock B, Gilbert J A, Tait K, et al. 2011. Bioturbation: Impact on the marine nitrogen cycle. Biochemical Society Transactions, 39(1): 315–320. doi: 10.1042/BST0390315
|
|
Li Bai, Cao Jie, Guan Lisha, et al. 2018. Estimating spatial non-stationary environmental effects on the distribution of species: a case study from American lobster in the Gulf of Maine. ICES Journal of Marine Science, 75(4): 1473–1482. doi: 10.1093/icesjms/fsy024
|
|
Li Mingkun, Zhang Chongliang, Xu Binduo, et al. 2020. A comparison of GAM and GWR in modelling spatial distribution of Japanese mantis shrimp (Oratosquilla oratoria) in coastal waters. Estuarine, Coastal and Shelf Science, 244: 106928
|
|
Liu Haiying, Wang Gui’e, Wang Xiuli. 2009. Genetic diversity analysis of mantis shrimp Oratosquilla oratoria from Dalian coast. Journal of Dalian Fisheries University (in Chinese), 24(4): 350–353
|
|
Lou Fangrui, Gao Tianxiang, Han Zhiqiang. 2019. Effect of salinity fluctuation on the transcriptome of the Japanese mantis shrimp Oratosquilla oratoria. International Journal of Biological Macromolecules, 140: 1202–1213. doi: 10.1016/j.ijbiomac.2019.08.223
|
|
Lu Binbin, Charlton M, Harris P, et al. 2014. Geographically weighted regression with a non-Euclidean distance metric: A case study using hedonic house price data. International Journal of Geographical Information Science, 28(4): 660–681. doi: 10.1080/13658816.2013.865739
|
|
Lui K K Y, Leung P T Y, Ng W C, et al. 2010. Genetic variation of Oratosquilla oratoria (Crustacea: Stomatopoda) across Hong Kong waters elucidated by mitochondrial DNA control region sequences. Journal of the Marine Biological Association of the United Kingdom, 90(3): 623–631. doi: 10.1017/S0025315409990841
|
|
Luo Chongxin, Lin Lei, Shi Jie, et al. 2021. Seasonal variations in the water residence time in the Bohai Sea using 3D hydrodynamic model study and the adjoint method. Ocean Dynamics, 71(2): 157–173. doi: 10.1007/s10236-020-01438-5
|
|
Manning R B. 1971. Keys to the species of Oratosquilla (Crustacea, Stomatopoda), with descriptions of two new species. Smithsonian Contributions to Zoology, 1: 1–16
|
|
Maynou F, Abelló P, Sartor P. 2004. A review of the fisheries biology of the mantis shrimp, Squilla mantis (L., 1758) (Stomatopoda, Squillidae) in the Mediterranean. Crustaceana, 77(9): 1081–1099. doi: 10.1163/1568540042900295
|
|
Nakajima M, Kodama K, Horiguchi T, et al. 2010. Impacts of shifts in spawning seasonality and size at maturation on the population growth of mantis shrimp in Tokyo Bay. Marine Ecology Progress Series, 418: 179–188. doi: 10.3354/meps08824
|
|
Ng Yingpei. 2013. The Ecology of Stomatopods in Matang waters with emphasis on Miyakea nepa and Oratosquillina perpensa [dissertation]. Kuala Lumpur: University of Malaya
|
|
Qiao Huijie, Feng Xiao, Escobar L E, et al. 2019. An evaluation of transferability of ecological niche models. Ecography, 42(3): 521–534. doi: 10.1111/ecog.03986
|
|
Ragheb E, Akel E S H K, Rizkalla S I. 2019. Analyses of the non-target catch from the Egyptian Mediterranean trawlers, off Port Said. The Egyptian Journal of Aquatic Research, 45(3): 239–246. doi: 10.1016/j.ejar.2019.07.003
|
|
Reaka M L. 1980. Geographic range, life history patterns, and body size in a guild of coral-dwelling mantis shrimps. Evolution, 34(5): 1019–1030
|
|
Reaka M L, Rodgers P J, Kudla A U. 2008. Patterns of biodiversity and endemism on Indo-West Pacific coral reefs. Proceedings of the National Academy of Sciences of the United States of America, 105(S1): 11474–11481
|
|
Robert K, Jones D O B, Roberts J M, et al. 2016. Improving predictive mapping of deep-water habitats: Considering multiple model outputs and ensemble techniques. Deep-Sea Research Part I: Oceanographic Research Papers, 113: 80–89. doi: 10.1016/j.dsr.2016.04.008
|
|
Tanaka K, Chen Y. 2015. Spatiotemporal variability of suitable habitat for american lobster (Homarus Americanus) in Long Island sound. Journal of Shellfish Research, 34(2): 531–543. doi: 10.2983/035.034.0238
|
|
Tao L S R, Lui K K Y, Lau E T C, et al. 2018. Trawl ban in a heavily exploited marine environment: Responses in population dynamics of four stomatopod species. Scientific Reports, 8(1): 17876. doi: 10.1038/s41598-018-35804-7
|
|
Wang Lei, Qiu Shengrao, Liu Shude, et al. 2020. Morphological difference analysis on three different stocks of Oratosquilla oratoria in the Bohai Sea and the Yellow Sea. Marine Fisheries (in Chinese), 42(6): 672–686
|
|
Wilson J. 2017. Learning, adaptation, and the complexity of human and natural interactions in the ocean. Ecology and Society, 22(2): 43. doi: 10.5751/ES-09356-220243
|
|
Wu Qiang, Guan Lisha, Shan Xiujuan, et al. 2019. Decadal variations in the community status of economically important invertebrates in the Bohai Sea. Acta Oceanologica Sinica, 38(10): 60–66. doi: 10.1007/s13131-019-1488-1
|
|
Xue Ying, Guan Lisha, Tanaka K, et al. 2017. Evaluating effects of rescaling and weighting data on habitat suitability modeling. Fisheries Research, 188: 84–94. doi: 10.1016/j.fishres.2016.12.001
|
|
Yang Jiming. 2001. A study on food and trophic levels of Bohai Sea Invertebrates. Modern Fisheries Information (in Chinese), 16(9): 8–16
|
|
Yang Mei, Li Xinzheng. 2018. Population genetic structure of the mantis shrimp Oratosquilla oratoria (Crustacea: Squillidae) in the Yellow Sea and East China Sea. Journal of Oceanology and Limnology, 36(3): 905–912. doi: 10.1007/s00343-018-7030-z
|
|
Yuan Ping, Wang Houjie, Wu Xiao, et al. 2020. Grain-size distribution of surface sediments in the Bohai Sea and the northern Yellow Sea: sediment supply and hydrodynamics. Journal of Ocean University of China, 19(3): 589–600. doi: 10.1007/s11802-020-4221-y
|
|
Zhang Yang, Han Zhiqiang, Gao Tianxiang, et al. 2018. Genetic structure analysis of mantis shrimp Oratosquilla oratoria based on mitochondrial DNA control region sequence. Genes and Genomics, 40(9): 1001–1009. doi: 10.1007/s13258-018-0707-z
|
|
Zhang Yunlei, Yu Huaming, Yu Haiqing, et al. 2020. Optimization of environmental variables in habitat suitability modeling for mantis shrimp Oratosquilla oratoria in the Haizhou Bay and adjacent waters. Acta Oceanologica Sinica, 39(6): 36–47. doi: 10.1007/s13131-020-1546-8
|
|
Zhao Wenju, Cao Taohong, Dou Pinxin, et al. 2019. Effect of various concentrations of superabsorbent polymers on soil particle-size distribution and evaporation with sand mulching. Scientific Reports, 9(1): 3511. doi: 10.1038/s41598-019-39412-x
|
Figures(7) / Tables(2)
Supported by:
Beijing Renhe Information Technology Co. Ltd
Lisha Guan, Xianshi Jin, Tao Yang, Xiujuan Shan. Ensemble habitat suitability modeling of stomatopods with Oratosquilla oratoria as an example[J]. Acta Oceanologica Sinica, 2023, 42(4): 93-102. doi: 10.1007/s13131-022-2051-z
| Grain size/mm | Grain size/μm | Wentworth class | Phi (Φ) range |
| 1 to 2 | 1000 to 2000 | very coarse sand | Φ=–1 to Φ=0 |
| 1/2 to 1 | 500 to 1000 | coarse sand | Φ=0 to Φ=1 |
| 1/4 to 1/2 | 250 to 500 | medium sand | Φ=1 to Φ=2 |
| 1/8 to 1/4 | 125 to 250 | fine sand | Φ=2 to Φ=3 |
| 1/16 to 1/8 | 62.5 to 125 | very fine sand | Φ=3 to Φ=4 |
| 1/32 to 1/16 | 31.25 to 62.5 | coarse silt | Φ=4 to Φ=5 |
| 1/64 to 1/32 | 15.625 to 31.25 | medium silt | Φ=5 to Φ=6 |
| 1/128 to 1/64 | 7.813 to 15.625 | fine silt | Φ=6 to Φ=7 |
| 1/256 to 1/128 | 3.906 to 7.813 | very fine silt | Φ=7 to Φ=8 |
| < 1/256 | < 3.906 | clay | Φ>8 |
DownLoad:
CSV
| Ord_Psand | Sand content | Ord_Pclay | Clay content | Ord_Mz | Mean size (Φ) | Ord_Sk | Skewness | Ord_Kt | Kurtosis |
| 1 | 0 to 10% | 1 | 0 to 5% | 1 | Φ<3 | 1 | –2.5 to –1.5 | 1 | 0.5 to 1.5 |
| 2 | 10% to 20% | 2 | 5% to 10% | 2 | Φ=3 to Φ=4 | 2 | –1.5 to –0.33 | 2 | 1.5 to 2.5 |
| 3 | 20% to 30% | 3 | 10% to 15% | 3 | Φ=4 to Φ=5 | 3 | –0.33 to 0.33 | 3 | 2.5 to 2.75 |
| 4 | 30% to 40% | 4 | 15% to 20% | 4 | Φ=5 to Φ=6 | 4 | 0.33 to 1.5 | 4 | 2.75 to 3 |
| 5 | 40% to 50% | 5 | 20% to 25% | 5 | Φ=6 to Φ=7 | 5 | 1.5 to 2.2 | 5 | 3 to 3.25 |
| 6 | 50% to 60% | 6 | 25% to 30% | 6 | Φ>7 | 6 | 2.2 to 3.5 | 6 | 3.25 to 3.5 |
| 7 | 60% to 70% | 7 | 30% to 35% | – | – | 7 | > 3.5 | 7 | 3.5 to 4.5 |
| 8 | 70% to 80% | 8 | 35% to 40% | – | – | – | – | 8 | 4.5 to 5.5 |
| − | − | − | − | − | − | − | − | 9 | >5.5 |
DownLoad:
CSV
| Grain size/mm | Grain size/μm | Wentworth class | Phi (Φ) range |
| 1 to 2 | 1000 to 2000 | very coarse sand | Φ=–1 to Φ=0 |
| 1/2 to 1 | 500 to 1000 | coarse sand | Φ=0 to Φ=1 |
| 1/4 to 1/2 | 250 to 500 | medium sand | Φ=1 to Φ=2 |
| 1/8 to 1/4 | 125 to 250 | fine sand | Φ=2 to Φ=3 |
| 1/16 to 1/8 | 62.5 to 125 | very fine sand | Φ=3 to Φ=4 |
| 1/32 to 1/16 | 31.25 to 62.5 | coarse silt | Φ=4 to Φ=5 |
| 1/64 to 1/32 | 15.625 to 31.25 | medium silt | Φ=5 to Φ=6 |
| 1/128 to 1/64 | 7.813 to 15.625 | fine silt | Φ=6 to Φ=7 |
| 1/256 to 1/128 | 3.906 to 7.813 | very fine silt | Φ=7 to Φ=8 |
| < 1/256 | < 3.906 | clay | Φ>8 |
| Ord_Psand | Sand content | Ord_Pclay | Clay content | Ord_Mz | Mean size (Φ) | Ord_Sk | Skewness | Ord_Kt | Kurtosis |
| 1 | 0 to 10% | 1 | 0 to 5% | 1 | Φ<3 | 1 | –2.5 to –1.5 | 1 | 0.5 to 1.5 |
| 2 | 10% to 20% | 2 | 5% to 10% | 2 | Φ=3 to Φ=4 | 2 | –1.5 to –0.33 | 2 | 1.5 to 2.5 |
| 3 | 20% to 30% | 3 | 10% to 15% | 3 | Φ=4 to Φ=5 | 3 | –0.33 to 0.33 | 3 | 2.5 to 2.75 |
| 4 | 30% to 40% | 4 | 15% to 20% | 4 | Φ=5 to Φ=6 | 4 | 0.33 to 1.5 | 4 | 2.75 to 3 |
| 5 | 40% to 50% | 5 | 20% to 25% | 5 | Φ=6 to Φ=7 | 5 | 1.5 to 2.2 | 5 | 3 to 3.25 |
| 6 | 50% to 60% | 6 | 25% to 30% | 6 | Φ>7 | 6 | 2.2 to 3.5 | 6 | 3.25 to 3.5 |
| 7 | 60% to 70% | 7 | 30% to 35% | – | – | 7 | > 3.5 | 7 | 3.5 to 4.5 |
| 8 | 70% to 80% | 8 | 35% to 40% | – | – | – | – | 8 | 4.5 to 5.5 |
| − | − | − | − | − | − | − | − | 9 | >5.5 |
DownLoad:
DownLoad:
