Optimization of environmental variables in habitat suitability modeling for mantis shrimp Oratosquilla oratoria in the Haizhou Bay and adjacent waters

Habitat refers to environmental conditions for species that provide essential living conditions for them (Shen et al., 2010). Once the suitable habitat shrinks or disappears, abundance of species is reduced, or species may even become extinct (Morrison et al., 2006; Xue et al., 2017). Ecosystem-based Fisheries Management (EBFM) has been widely considered as an essential approach for addressing the crisis of fishery management around the world. However, EBFM is generally hampered by a lack of information about habitats and other ecosystem components (Pikitch et al., 2004). Understanding the spatiotemporal dynamics of species distributions and its relationships with environmental conditions are crucial for the EBFM and the evaluation of the effects of climate change (Liu et al., 2019). Therefore, the unerring analysis and evaluation of habitat suitability of species are important prerequisites for EBFM (Gong et al., 2011; Tanaka and Chen, 2016). Habitat suitability index (HSI) models have been widely used in analyzing the relationship between species abundance and environmental factors, and ultimately informing management of wildlife (Chang et al., 2012; Tanaka and Chen, 2015). HSI models have recently been applied in the assessment and forecast of suitable habitat for fishery resources (Tian et al., 2009; Chen et al., 2010; Tanaka and Chen, 2015; Zou et al., 2016). The suitability indices (SIs) are based on functional relationships between several key environmental variables and species abundance to evaluate the habitat suitability of target species (Brown et al., 2000; Tian et al., 2009; Chen et al., 2010). All of the SIs are combined to a composite HSI score, ranging from 0 (“poor”) to 1 (“good”) to indicate the habitat quality (Brooks, 1997).

Typically, all environmental variables are modelled with equal weights in HSIs (Vayghan et al., 2013; Yu et al., 2016). However, the response of species abundance to each environmental variable is different, and habitat requirements may change over life history stages and seasons (Valavanis et al., 2004). Therefore, it is necessary to select the optimal combination of environmental variables to build HSI models. Many statistical models are widely used to analyze the relationships between species abundance and habitat suitability (Ahmadi-Nedushan et al., 2006; Xue et al., 2018). The statistical non-linear, generalized additive model (GAM), is one of the most suitable models to explore the relationships between fish distribution and environmental factors (Hastie and Tibshirani, 1990; Leathwick et al., 2006; Ptacnik et al., 2008; Chang et al., 2010; Schmiing et al., 2013). The advantage of GAM is that it offers a flexible and robust approach to select environmental variables to include in HSI models. GAMs were used to determine which environmental variables to be included in the HSI models. This may be attributed to the flexibility of GAMs with their semi-parametic functions being able to accommodate non-linear relationships between abundance and environmental variables (Becker et al., 2010).

The relative importance of each environmental variable can be determined by its contribution to the total variance explained by the BRT model (Xue et al., 2017). Boosted regression trees (BRTs) use an integrated learning method based on decision tree theory (Hastie et al., 2001; Gao et al., 2015), which has been applied in studies of fish habitat (Compton et al., 2012; Lewin et al., 2014; Xue et al., 2017). BRT models are based on two statistical algorithms (i.e., regression trees and boosting), and incorporate the advantages of both methods, which improve the predictive performance of habitat modeling (Elith et al., 2008). The regression trees are capable of handling missing values and can automatically fit the interactions between predictors (Compton et al., 2012). Boosting improves model performance by focusing on the observations that are difficult to predict with a sequential model fitting process and by adding a probabilistic component to optimize predictive performance (Compton et al., 2012). BRTs are able to deal with the non-linear, correlated and interacting variables (Torres et al., 2015). The BRT model is effective in determining the relative importance of environmental variables in HSI modeling (Xue et al., 2017).

The Haizhou Bay is a trumpet-shaped open coastal bay lying on the western margin of the South Yellow Sea. This bay is rich in commercially important fishes and other species (Su et al., 2015). It covers an area of approximately 876.39 km² (Zhang et al., 2014). The length of the coastline is 86.81 km, and the bay’s maximum width is approximately 42 km. In recent years, a tremendous decline in fishery resources and diversity has taken place in the Haizhou Bay as a result of over-fishing, coastal development and environmental pollution (Chen, 1993; Tang et al., 2011; Su et al., 2015). Mantis shrimp (Oratosquilla oratoria) is now one of the dominant species in this area (Wang et al., 1996b; Xu et al., 2017). It can be found in a cave of sediment and lives in shallow coastal water (<60 m), distributed from the Gulf of Pebbly in Russia to the coasts of Japan, China, Philippines, Malay Peninsula and the Hawaiian Islands (Wang et al., 1998). In particular, mantis shrimp occurs in a wide range of continuously distributed habitats, extending broadly from the inshore areas down to coral reefs along the China coastal waters (Du, 1993; Huang, 2008). However, there is a lack of information available on the habitat suitability of mantis shrimp in China seas.

In this study, HSI models of mantis shrimp (Oratosquilla oratoria) were established based on the bottom trawl surveys in the Haizhou Bay and adjacent waters during spring and fall in 2011 and 2013–2017. The goals of this study were to determine the optimal combination of environmental factors in HSI modeling based on GAMs and BRTs, and examine the seasonal and ontogenetic variations in the suitable habitat of mantis shrimp. Our study provides a better approach for the evaluation of suitable habitats for marine organisms and a scientific basis for the protection and sustainable utilization of mantis shrimp resources.

Mantis shrimp samples were collected from bottom trawl surveys in the Haizhou Bay and adjacent areas during spring (May) and fall (September to October) in 2011 and 2013–2017 (no survey in 2012) (Fig. 1). The survey area was 34°20′–35°40′N, 119°20′–121°10′E. Surveys were designed by stratified random sampling (Xu et al., 2015). According to the differences in oceanographic, geological and biological characteristics of the survey area, a certain number of stations were randomly selected in each stratification. A target of 24 stations were selected in 2011, and the sampling stations were optimized to 18 after 2013 (Xu et al., 2015). The bottom trawl was towed at a speed of 2–3 kn for aiming for 1 h, with a trawl net of 12 m width and mesh size of 17 mm. There were two bottom trawl vessels with the same size, power and trawl net in our study, and tows were only conducted during daytime. Further details about the survey design, gears and optimization process are available in Xu et al. (2015). At each survey station, a CTD system (XR-420) was used to measure environmental data including depth, bottom temperature and bottom salinity. Distance offshore was the distance between sampling station and the nearest coastline. The sediment data of the Haizhou Bay were provided by the College of Environmental Science and Engineering at Ocean University of China (Li et al., 2014). Ten and eight kinds of sediment types were detected in spring (CS: coarse sand; MCS: middle coarse sand; S: sand; MFS: middle fine sand; TS: silty sand; YS: clayey sand; ST: sandy silt; MT: muddy silt; STY: sand-silt-clay; and TY: silty clay) and fall (CS, S, MFS, TS, YS, ST, STY and TY), respectively.

Figure 1.  Study area and sampled stations (dots) by bottom trawl surveys in the Haizhou Bay of 2011 and 2013–2016 (all years combined).

Before data analysis, the survey data were standardized by 2 kn and 1 h, and the catch per unit area (ind./km²) was used as the abundance index. The survey data were analyzed separately by seasons (spring and fall) and ontogenetic stages (juvenile: ≤80 mm body length, adult: >80 mm body length). Length at 50% maturity (L50%) for mantis shrimp is 80 mm (Hamano and Matsuura, 1984; Ohtiomi et al., 1988; Xu et al., 1996).

The environmental data (depth, bottom temperature and bottom salinity) extracted from the FVCOM model was used to predict HSI of mantis shrimp in spring and fall of 2011 and 2013–2017. The Finite-Volume Community Ocean Model (FVCOM) is a model of the coastal ocean circulation developed jointly by the University of Massachusetts Dartmouth (UMASSD) and the Woods Hole Oceanographic Institution (WHOI), which is an unstructured-grid, finite-volume, free-surface, 3-D primitive model (http://fvcom.smast.umassd.edu/). The model has been widely used for coastal environment predictions and management (Chen et al., 2016; Li et al., 2017).

Previous studies showed that several environmental factors, including bottom temperature, bottom salinity, depth, distance offshore and sediment type, were important for the distribution of mantis shrimp (Liu et al., 2014). Prior to selecting the environmental variables, all candidate variables were tested by variance inflation factors (VIF) in R statistical program, with $\sqrt {\rm{VIF}} > 2$ being an indication of strong collinearity and should be dropped in the model (Kabacoff, 2015). The VIF value is determined using vif function in “car” package of R software (version: 3.4.2) (Fox and Weisberg, 2019).

GAMs were used to determine which environmental variables to be included in the HSI models. GAM is a generalization of the usual linear regression model and additive in the predictor effects. GAM is usually expressed as:

where $ {f_i}\left( {{x_i}} \right)$ is a smoothing spline function of explanatory variables, α is the intercept term, n is the number of predictors, and ε is the residual error term (Hastie and Tibshirani, 1990), which gives GAMs a partially non-parametric aspect. The log-transformed mantis shrimp abundance lg(Y) was used as the response variable with Gaussian distribution. In order to avoid overfitting, all the smoothers were constrained with 4 kn. The stepwise variable selection procedure was used to select significant variables, which started with a null model and added one predictive variable to the present model at each step (Luan et al., 2018). Akaike Information Criterion (AIC) (Akaike, 1998) was used in variable selection for GAM, in which variables in GAM with the lowest AIC value was regarded as the significant variables.

BRTs were used to determine the weight of each environmental variable in HSI models based on their relative contributions (Chang et al., 2010; Yi et al., 2016). The boosting algorithm is optimized by learning rate and tree complexity (Elith et al., 2008). The learning rate calculates the residual variation in tree construction, and tree complexity estimates the interaction between the predictor variables (Compton et al., 2012). In our study, tree complexity and learning rate were 5 and 0.01, respectively. The relative importance of a predictor in the BRT model was determined by its contribution to the model, measured by the number of times the tree was split by that environmental variable. Each environmental variable was assigned a weight value by its relative contribution (%) to the total deviation explained by the BRT model. The BRT model was performed using the “gbm.step” function of the “gbm” package in an R programming environment (Ridgeway, 2015).

HSI models for each season and ontogenetic stage of mantis shrimp were built. The spline smooth regression was used to fit the relationship between each environmental variable and suitability index (SI) (Chang et al., 2012). The value of SI ranges from 0 to 1, and the formula is as follows:

where $ \hat Y$ is the predicted abundance in spline smooth regression, and $ {\hat Y_{{\rm{max}}}}$ and $ {\hat Y_{{\rm{min}}}}$ are the maximum and minimum of the predicted abundance, respectively. The suitability index is the highest with a SI value of 1, while the lowest abundance has the lowest suitability index of 0. SI values between 0.7 and 1.0 are considered suitable habitat range (Brooks, 1997; Xue et al., 2017).

HSI is a comprehensive model composed of SI values. The value of HSI is between 0 and 1, indicating “poor” to “good” habitat quality. In this study, we used two empirical formulas, arithmetic mean model (AMM) and geometric mean model (GMM), to build HSI models, and assigned the weight based on BRT output to each environmental variable. The formulas are as follows:

(1) Weighted arithmetic mean model (AMM) (Wakeley, 1988):

(2) Weighted geometric mean model (GMM) (Vincenzi et al., 2007):

where SI_i is the SI value of environmental variable i, w_i is the weight assigned to environmental variable i based on BRT model output, and n is the number of environmental variables.

HSI models were constructed based on observed environmental data from 2011 and 2013–2016. Cross-validation was applied to assess the performance of HSI models. A total of 80% of the data subset (training data) were randomly selected for model development, and the remaining 20% of the data (testing data) were used to evaluate the performance of the model (Smith, 1994; Zuur et al., 2007). The regression tests and Akaike information criterion (AIC) were also used to evaluate the performance of HSI models. The model that yielded the minimum Akaike information criterion (AIC) value and maximum R² was selected as the most suitable model. Four types of HSI models were evaluated in our study, i.e., (1) the non-optimized models, (2) the BRT informed HSI model (only optimized with weighted environmental variables), (3) the GAM informed HSI model (only optimized with selected environmental variables), and (4) both BRT and GAM informed HSI models. A total of 100 regression tests were performed for each model. The predicted environmental data extracted from the FVCOM model were used to make the maps of HSI values of mantis shrimp. The predicted HSI distribution in 2017 was then overlaid and tested with the observed abundance of mantis shrimp (log-transformed) in 2017 by Pearson correlation tests to verify whether there was a high consistency between them (Xue et al., 2017; Yu et al., 2018, 2019). HSI maps were constructed using ordinary kriging (Xue et al., 2017).

The degree of multicollinearity was measured by VIF, with $ \sqrt {\rm{VIF}} $ values of candidate variables of the four groups being all less than 2. Therefore, all factors were retained for both GAMs and BRT models.

AIC values in GAMs with the combination of depth, bottom temperature, distance offshore and sediment type were the smallest for all groups (spring-juvenile, spring-adult, fall-juvenile and fall-adult) (Table 1). Bottom salinity was not selected by the GAMs. In addition, the cumulative deviance explained in four groups were 58.96%, 62.22%, 51.90% and 53.43%, respectively (Table 1).

In spring, the contribution of depth to the total variation was the highest for juveniles (43.70%) and adults (66.18%), followed by bottom temperature and sediment type. The contribution of distance offshore was lowest, at 1.72% and 1.40% for juveniles and adults, respectively (Fig. 2).

Figure 2.  Relative contribution of different environmental variables to the total deviance explained by the boosted regression tree (BRT) models for different groups of mantis shrimp (Oratosquilla oratoria) in the Haizhou Bay.

In fall, the contribution of depth to the total variation was the highest, accounting for 73.95% and 49.14% of total variance for juveniles and adults, respectively. The factor with the second highest contribution was sediment type for juveniles (11.56%) and bottom temperature for adults (16.19%). The contribution of distance offshore to the total deviation explained was the lowest for both life history stages, being 4.30% for juveniles and 13.70% for adults (Fig. 2).

The optimal ranges of environmental factors of mantis shrimp (SI>0.7) varied by seasons and life history stage (Fig. 3). The optimum range of bottom temperature was from 16.0°C to 18.0°C for both juveniles and adults in spring, while the optimum ranges of bottom temperature for juveniles and adults in fall were from 17.8°C to 18.3°C and 17.8°C to 21.2°C, respectively. The suitable depth was the same in both seasons, with juveniles and adults found in shallow coastal waters less than 10 m in depth. Suitability indices decreased gradually with depth. The optimal range of distance offshore was less than 10 km for all groups. Mantis shrimp preferred to inhabit sediment types of silty clay and clayey sand in spring, and silty sand in fall (Fig. 3).

Figure 3.  Suitability index (SI) curves of bottom temperature, depth, distance offshore and sediment types for spring-juvenile (a), spring-adult (b), fall-juvenile (c) and fall-adult (d) of mantis shrimp (Oratosquilla oratoria). The dashed line (SI=0.7) indicated the optimal range of environmental variables. CS: coarse sand; MCS: middle coarse sand; S: sand; MFS: middle fine sand; TS: silty sand; YS: clayey sand; ST: sandy silt; MT: muddy silt; STY: sand-silt-clay; TY: silty clay. Sediment types were ordered by sediment grain size (from big to small).

AMM HSI models informed by both BRT and GAM performed better in the cross-validation test compared with the other three models (non-optimized model, BRT informed HSI model and GAM informed HSI model), with the highest mean R² and the lowest AIC values for all groups (spring-juvenile, spring-adult, fall-juvenile and fall-adult) (Table 2). In addition, BRT informed HSI model showed better performance than GAM informed HSI model (Table 2).

The predicted HSI values in 2017 were overlaid with the observed abundance of mantis shrimp in 2017 to further evaluate the performance of the HSI models (Fig. 4). In spring, the non-optimized and optimized HSI models (both BRT and GAM informed HSI model) predicted similar distribution patterns of habitat suitability for all groups, with higher HSI values mainly distributed in the coastal areas of the Haizhou Bay (at depths less than 20 m). However, in the fall, the optimized HSI matched better with the actual catches of all groups than those derived from the non-optimized HSI models (Fig. 4), which are consistent with the results of correlation test (Table 3). In addition, the maximum HSI values for non-optimized (Figs 4a-d) and optimized (Figs 4e-h) HSI models were about 0.65 and 1.00, respectively. The non-optimized HSI models tend to underestimate HSI values and the range of optimal habitats (Fig. 4).

Figure 4.  Spatial distribution of observed abundance (black circles) of mantis shrimp in 2017 overlaid with the predicted HSI values in 2017 (colour contours) based on non-optimized (a–d) and optimized (e–h) HSI models.

The suitable habitats of mantis shrimp in the Haizhou Bay and adjacent waters was mainly distributed in shallow waters less than 20 m and differed by seasons (Fig. 5). There were no obvious ontogenetic variations. In spring, the habitat suitability in areas south of 35°N was higher than that in the north. The areas deeper than 20 m were less suitable, with HSI values being lower than 0.4. In fall, mantis shrimp were distributed throughout the survey area, and the HSI values in depths more than 20 m were significantly higher (P<0.05) than those in spring. In addition, there were some oval patches of lower HSI in fall (Fig. 5).

Figure 5.  Spatial distribution of the mean habitat suitability index (HSI) values for spring-juvenile (a), spring-adult (b), fall-juvenile (c) and fall-adult (d) groups of mantis shrimp based on optimized HSI model in the Haizhou Bay and adjacent waters during 2011 and 2013–2016.

In the past, all relevant environmental variables with the same weight were incorporated in HSI modeling due to the lack of information on the relative importance of each environmental variable for a species (Vayghan et al., 2013; Yu et al., 2016). However, environmental variables affecting distribution of species are usually different in different areas, seasons and even different growth stages for species (Valavanis et al., 2004). Therefore, it is necessary to select the best combination of environmental variables to be included in HSI modeling. In this study, we developed a framework to deal with this issue. The results (Table 2) showed that both BRT and GAM informed HSI models performed better than the other three HSI models, and BRT informed HSI model showed better performance than GAM informed HSI model. For optimization of environmental variables in HSI modeling, weighting was more important than selection.

Understanding the response of species to environmental factors furthers our understanding on physiological and behavioral characteristics of different species (Zeng and Yeo, 2018). For selection of candidate environmental factors in modeling, Chang et al. (2010) found that bottom temperature, salinity, longitude, latitude, depth and distance offshore were important environmental variables affecting the spatial distribution of American lobster (Homarus americanus). In addition, there were differences in environmental variables for different seasons, sex and life history stages of this species. The selection of environmental variables is an important factor affecting the performance of HSI models (Chen et al., 2008; Tian et al., 2009; Xue et al., 2017).

Marine organisms respond differently to a variety of ecological factors (Valavanis et al., 2004). In this study, five environmental variables including bottom temperature, bottom salinity, depth, distance offshore and sediment type, were chosen as potential factors for HSI modeling of mantis shrimp. Following the GAM results, bottom salinity was excluded from the HSI models for all groups of mantis shrimp. Previous studies suggested that mantis shrimp can live in a wide range of bottom salinity (Wang et al., 1996b), with suitable salinity between 24 and 36 (Liu et al., 2006). Based on our survey data, bottom salinity in the Haizhou Bay ranged from 28.36 to 32.03 in spring and from 27.30 to 31.94 in fall during 2011, and 2013–2017. Therefore, the salinity was suitable for mantis shrimp in this area and it was not a main factor that restricted mantis shrimp distribution in the Haizhou Bay.

Different environmental variables usually do not have the same impacts in HSI modeling (Xue et al., 2017; Yu et al., 2018). Therefore, weight for each factor should be considered in HSI modeling. In previous studies, the weight of each environment variable was determined by expert judgment or mathematical statistics (Chen et al., 2011; Liang et al., 2015; Yi et al., 2016; Xue et al., 2017). For example, Gong et al. (2012) proposed ten different weights for different environmental variables for neon flying squid (Ommastrephes bartramii) in the Northwest Pacific and evaluated the impact of different weights on HSI models. In particular, SST (sea surface temperature, 0.33–0.7) was given the highest weight and the weights for GSST (horizontal gradient of SST, 0.33–0.15) and SSH (sea surface height, 0.33–0.15) varied with the weight for SST. Liang et al. (2015) assigned the weight of the environmental variables based on statistics method, canonical correspondence analysis (CCA), to assess the suitability of waterfowl habitats in the eastern Dongting Lake. In this study, the weights of environmental factors were determined by the BRT model and varied with seasons (spring and fall) and ontogenetic stages (juvenile and adult). The performance of weighted HSI models was better than the unweighted model. Since BRTs are able to accommodate different types of predictors and missing values, automatically fitting interactions between predictors (Friedman and Meulman, 2003), the performance of HSI models was improved through application of appropriate weights to environmental variables. In this study, GAMs selected the same variable, and BRTs allocated different weights for four modeling environmental variables in different seasons and growth stages, which is related to the ecological habits of the species and their adaptability to the environment.

Wu et al. (2015) found that the abundance of mantis shrimp was the highest in the coastal shallow waters of China in May. Mantis shrimp spawns in shallow waters in spring and migrates short distances to deep waters in fall (Wu et al., 2015). The habitat suitability results from this study were similar with those from previous studies (Wang et al., 1996a). Depth was the most important factor in all mantis shrimp groups. Depth and other spatial variables such as longitude and latitude, practically serve as proxies for other important variables (Xue et al., 2017). Distance offshore was used as spatial variable instead of longitude and latitude in our study. Depth can affect various other factors, such as transparency (most often refer to transparency and translucency, the physical property of allowing the transmission of light through a material), water flow, dissolved oxygen, and food availability (Chen, 2004).

Temperature and sediment are important ecological factors that regulate the growth, development and reproduction of many marine organisms (Shen, 2010). Mantis shrimp are mainly distributed in coastal areas with high temperature (14.0–20.7°C; Liu et al., 2014). The spawning season of mantis shrimp is between May and July (Deng et al., 1988; Wang et al., 1998). Compared to other months, temperature from May to July is more suitable for spawning adults (Wang et al., 1998). After the spawning season, there was no significant difference in suitability indices of temperature between juveniles and adults (P>0.05). It prefers to live in the silt and sand habitats (Wang et al., 1996b). In this study, suitable habitats for mantis shrimp included sediment of TS and TY probably because mantis shrimp can easily build caves and spawn in these substrates (Wang et al., 1998; Wu et al., 2015). Therefore, substrate types should be included in HSI modeling for mantis shrimp.

In spring, the Yellow Sea Cold Water Mass in Shandong Peninsula restricts the distribution of mantis shrimp in coastal areas, and the cold water zone in the Yellow Sea is a limiting factor for the distribution of mantis shrimp within the 20 m isobath (Zhao et al., 1992; Su and Huang, 1995; Shi, 2019). As the temperature increases in shallow waters of the Haizhou Bay, the nutrients and prey organisms are sufficient (Pang et al., 2015), which provides enough food for the growth of mantis shrimp. After September, the Yellow Sea Cold Water Mass gradually becomes weaker, and the distribution area of mantis shrimp increases (Shi, 2019). Mantis shrimp move to deeper waters for over-wintering in fall when the temperature in shallow coastal waters begins to decline (Tang and Ye, 1990). Therefore, there are seasonal variations in the suitable habitats of mantis shrimp in the Haizhou Bay.

Among four variables in modeling (bottom temperature, depth, distance offshore, and sediment type), only bottom temperature changed over the year. In spring, due to samples were collected in May, the bottom temperature remained basically stable. Therefore, there was no obvious interannual variation for HSI values in spring. In fall, the change of bottom temperature between September and October was slightly larger than that in spring, so the average HSI value in fall from 2011 to 2016 (except 2012) was different from the HSI map in 2017 in some areas. Interestingly, there were some oval patches of lower HSI in map of mean HSI in fall. It also could be attributed to slight change in the bottom temperature in different sampling months among these years in fall.

Fishery managers hope to obtain a reliable habitat distribution map to determine the key habitats of species, which is important for the protection of target species (Vinagre et al., 2006; Chang et al., 2012). Information on the spatial distribution of fishery organisms is important in the management and monitoring of fishery resources. It is necessary to examine the spatial distribution of key marine species and their habitats for the Ecosystem-based Fishery Management (EBFM) and the evaluation of the effects of climate change (Liu et al., 2019). HSI models are important methods in understanding species distribution and analyzing the relationship between species abundance and environmental factors (Johnson et al., 2013). The development of marine protected areas (MPAs) and EBFM requires to unerringly asses the habitat of target species and their relationship with environmental factors, especially in areas under severe climate changes and environmental pollution (Xue et al., 2017).

In this study, the non-optimized HSI models may lead to underestimation of HSI values and the range of optimal habitats, resulting in a conservative estimation of optimal habitats and potential MPAs. Optimized HSI models can provide more reliable estimates of suitable habitats for marine organisms and provide a better method for habitat research. This approach could be extended to other marine organisms to enhance understanding of the habitat suitability of target species. Further studies are needed to improve the accuracy of model prediction and explore the sensitivity of the predicted HSI values to the FVCOM output. In addition, more biotic and abiotic predictors (such as predators, competitors, oxygen, fishing pressure and hydrodynamic forces) should be considered in the evaluation of the model performance (Planque et al., 2011; Cormon et al., 2014).

We are grateful to many colleagues and graduate students in the laboratory of Fisheries Ecosystem Monitoring and Assessment for their assistance in field sampling and sample analyses. We thank three anonymous reviewers for their detailed and constructive comments which greatly improve this paper.