Assessing the Distribution and Sustainable Exploitation of Lophius litulon in Marine Areas Off Shandong, China

In recent years, the proportion and economic value of Lophius litulon (family Lophiidae) in the coastal fishery off Shandong Province, China has increased. In this study, we mapped the distribution of L. litulon abundance [catch per unit effort (CPUE)] and applied a generalized additive model (GAM) to explore the relationship between CPUE and environmental factors. Two data-limited methods (the BSM related to the CMSY method and the AMSY method) were used to evaluate the stock status and relevant fishery reference points of L. litulon. The results showed that the L. litulon stock was concentrated in the central Yellow Sea, at 34.0°–37.0° N and 121.0°–124.0° E, and the highest average CPUE of L. litulon in this area occurred in winter. The three most significant environmental factors affecting species abundance were bottom temperature, bottom salinity, and depth. L. litulon was most abundant when bottom temperature ranged from 5.8 to 10.6°C, depth was > 18 m, and bottom salinity varied from 31.0 to 33.2‰. The BSM and AMSY models indicated that the L. litulon stock was unhealthy and had been overfished in recent years, as its biomass remained below the level that can support maximum sustainable yield. The relative exploitation ratios were also high. These results provide the basic data for improving sustainability of the exploitation of L. litulon in the Yellow and Bohai Seas.


INTRODUCTION
Marine fishery resources face a variety of threats globally, which has roused international concern (Ricard et al., 2012). Multiple stresses, such as overfishing, global warming, and pollution have changed the structure of fishery resources and caused a decline in catch quantity and quality (Watson et al., 2013;Costello et al., 2016). Currently, approximately 34.2% of fish stocks are caught at unsustainable levels (FAO, 2020). China has been one of the largest national fisheries in the world since 1989 and has witnessed a decline in its coastal fishery resources (Wang et al., 2006;Zhang and Qiu, 2019). Shandong Province, adjacent to the Yellow Sea and Bohai Sea, is one of the major fishing provinces in China. Shandong coastal waters are the most productive areas for fishing, as they are spawning and feeding grounds for a variety of commercial fish. The fishing boats in Shandong coastal waters are mainly from Shandong province, but also from Liaoning, Jiangsu, Hebei Province, and Tianjin city. According to the 2010-2020 China Fishery Statistical Yearbook, the number of domestic offshore fishing vessels in Shandong Province in 2020 decreased by 8.4% compared with 2019, and by 23.8% compared with 2010, indicating a continuous decline (Fisheries and Fisheries Administration, 2010, 2020. The marine fishery yield of Shandong Province mainly comes from the Yellow Sea and Bohai Sea, which provides more than 80% of the total catch, with the remainder mainly from the East China Sea. In the past few decades, Shandong Province has benefited substantially from marine fishery development, at the cost of resources and environmental deterioration (Chen, 1991). Pollution in coastal waters, excessive fishing intensity, and a decline in food resources have reduced catches of high value species, such as Trichiurus lepturus (Zou et al., 2019). Species with traditionally less economic value comprise as much as 60%-70% of the current catch, thus lowering the overall fishing benefit (Huang, 2012). In this context, changes in the distribution and resource status of marine organisms require urgent study to ensure current exploitation are sustainable.
Lophius litulon (family Lophiidae) is a demersal marine fish that is distributed in the northwestern Pacific constituting an important part of the catches of neighboring countries. Ji et al. (2007) showed that there was no significant difference in genetic structure between two geographical populations of L. litulon in the Yellow Sea and the Sea of Japan. Li X. et al. (2021) showed that L. litulon was the dominant species among 134 captured fish in the Shandong coastal waters. Furthermore, in waters off southwestern Korea, L. litulon was the dominant species, accounting for more than 60% of the total catch, together with Pampus echinogaster, Trichiurus lepturus, Engraulis japonicus, and Larimichthys polyactis (Kim et al., 2007). With the decline of traditional commercial fish resources, L. litulon, whose economic value has increased in recent years, has become a common species in China's domestic fish market and one of the main fish products exported abroad in recent years. Most previous studies on L. litulon have focused on its biological characteristics (Xu et al., 2010;Zhang et al., 2011), feeding habits (Xue et al., 2007), and migratory distribution (Michiol et al., 2002). Assessments of resource status and the influence of environmental factors are scarce. Although China has adopted a series of policies to protect coastal fisheries, the L. litulon population in the Yellow Sea is thought to have a simpler age structure and smaller body size (Li et al., 2012). Sustainable exploitation of L. litulon is required.
The assessment of fishery resources can be used as a scientific basis for fishery management, but few assessments of the L. litulon stock in the Yellow Sea have been conducted due to the lack of data (Wang et al., 2020). In recent years, three computerintensive methods, the Monte Carlo method CMSY (catchmaximum sustainable yield), the related Bayesian Schaefer model (BSM), and the abundance maximum sustainable yield (AMSY) have been proposed to evaluate stocks and related reference points of fishery resources in data-poor situations (Froese et al., 2017(Froese et al., , 2020. In particular, CMSY uses only a time series of catches and ancillary qualitative information to quantify the stock status and related fisheries reference points. In cases where relative abundance data [i.e., catch per unit effort; CPUE] are available in addition to catch data, the BSM method can be used to combine information from both datasets. AMSY is the most recent method for assessing fish populations based on abundance (CPUE) time series (Froese et al., 2020). Froese et al. (2018) and Palomares et al. (2018) showed that the fish stock status could be defined based on the B/B MSY and F/F MSY in the final year of the time series. When B/B MSY ≥ 1 and F/F MSY ≤ 1, the assessed stock is in a healthy state; when B/B MSY < 1 and F/F MSY > 1, the stock is overfished, with higher F/F MSY and lower B/B MSY indicating more severe overfishing.
A complex functional relationship exists between CPUE and related influential factors. This relationship can be modeled using species distribution models, such as the generalized additive model (GAM) (Hastie and Tibshirani, 1990;Xiao et al., 2004). In recent years, the GAM has been widely used to study the association between the temporal and spatial distribution of fishery resources, environmental factors, and fish stocks Ma et al., 2021). For example, Ma et al. (2021) used GAM to analyze the relationship between nominal and standardized CPUE and environmental factors for L. litulon in the Yellow Sea and Bohai Sea. Li et al. (2012) showed that the relationship between the spatial and temporal distribution of CPUE and environmental factors (such as year, position, water depth, and sea surface temperature) for L. litulon in the southern Yellow Sea was better explained by a GAM than a generalized linear model.
This study illustrated the temporal and spatial distribution of L. litulon CPUE in the Yellow Sea and Bohai Sea based on the production statistics from fishing vessels in Shandong Province from 2012 to 2019 and bottom trawl survey data in Shandong Province from 2016 to 2017. We evaluated the stock status and related fishery reference points using BSM and AMSY methods, and analyzed the relationship between CPUE and marine environmental factors such as surface temperature, salinity, and depth. The results can inform the Maritime Shandong Strategy put forward by Shandong Province and help fisheries advance into the era of farming the sea, herding, and fishing. These methods also provide theoretical baseline information for the management and development of L. litulon resources in the Yellow Sea and Bohai Sea and promote the sustainable exploitation of marine fishery resources.

Study Area
The Shandong Peninsula protrudes into the Bohai Sea and Yellow Sea, with numerous rivers flowing into both. The coastal waters of Shandong provide breeding, feeding, and nursery habitats for many fishery resources in the Yellow Sea and Bohai Sea, and support rich fishery resources for neighboring countries . The distribution of fishing areas for L. litulon in this study is shown in Figure 1.

Datasets
The L. litulon datasets used for the BSM and AMSY models were derived from the fishing logs of offshore fishing vessels in Shandong Province from 2012 to 2019, which included fishing power, mode of operation, fishing area, operating time, catch species, and yield. However, no records of fishing areas were available in the logbooks for 2017. Therefore, after standardizing the annual mean nominal CPUE for the other years, interpolation was used to obtain the mean nominal CPUE for 2017. The spatial resolution of each fishing area was 0.5 • × 0.5 • , and its location was represented by the latitude and longitude of the center point.
Commercial fishery monitoring and trawl surveys were combined to obtain the catch rate of L. litulon for each fishing area (Pecquerie et al., 2004;Gonzalez et al., 2021). The datasets used to draw monthly spatial distribution maps of L. litulon CPUE and to fit the GAM were obtained from (i) fishing logs of commercial fleets from 2014 to 2016, and (ii) bottom trawl survey data for the Shandong inshore fishery resources in October 2016 and January, May, and August 2017. The fishing data of L. litulon from 2014 to 2016 were used, as they were complete and more representative. The monthly mean nominal CPUE was calculated using both sources of data. Mean values represented nine months of the year, namely January-May and September-December, as data were not available from June to August due to the closed fishing season.
Environmental data from dataset "ii" (the bottom trawl resource survey) were used to fit the relationship between L. litulon and environmental factors. Environmental factors used in this study included surface temperature (ST), bottom temperature (BT), surface salinity (SS), bottom salinity (BS), water depth (D), and surface chlorophyll a concentration (SChl_A). Missing SChl_A data were supplemented by the MODIS_Aqua model on the NASA Ocean Color website 1 , and the SS and ST data were supplemented by the AQUA_MODIS model on the Physical Oceanography Distributed Active Archive Center website 2 . The spatial resolution of the environmental data was higher than that of the fishery data, so this study used the average value of the environmental data in each fishing area to match the spatial-temporal resolution of the L. litulon fishery datasets.

Calculation of Nominal Catch Per Unit Effort
In this study, the nominal CPUE [CPUE j , kg/(kW * d)] was calculated as follows (Maunder and Langley, 2004): where, C ij (kg) is the total monthly catch of fishing vessel i in fishing area j, P i (kW) is the power of operating fishing vessel i, d ij is the number of days of fishing vessel i operating in fishing area j in a month, and n j is the number of all fishing vessels operating in fishing area j in a month. The resulting nominal CPUE (and catch) datasets for commercial vessel logs and trawl survey are detailed in Supplementary Tables 1, 2.

Distribution Maps
Spatial distribution contour maps of L. litulon abundance were drawn by the ordinary kriging method using Surfer16 (Golden Software, Colorado). The data used to map the species distribution are presented in Supplementary Tables 2, 3, and the annual/monthly mean nominal CPUE is given in Supplementary Tables 4, 5.

Generalized Additive Model Method
A generalized additive model (GAM) was used to express the nonlinear relationship between the relative abundance of L. litulon and various environmental factors. GAM is presented as follows (Hastie and Tibshirani, 1990;R Core Team, 2018): where, g() is the connection function, Y is the abundance of L. litulon, β is the intercept term (as the response variable), f n (x n ) is the nonparametric function used to describe the relationship between g(Y) and the nth explanatory variable, which is estimated by the spline smoothing function, n is the number of selected environmental variables, and ε is the random error term (Xiao et al., 2004).
The connection function varies according to the actual distribution of response variable Y. In this study, the connection function g(Y) = log(CPUE+1) was used as the response variable, and month, ST, BT, SS, BS, D, SChl_A, Lon, and Lat were used as explanatory variables. Month was classified as a discrete variable, and all other variables were classed as continuous. The error distribution of the model was assumed to have a Gaussian distribution. Model construction was conducted using the RStudio software.

Bayesian Schaefer Model Method
In addition to the time series of catch and abundance data, BSM also requires ancillary qualitative information, that is, priors for relative biomass, intrinsic rate of population increase (r), and carrying capacity (k). According to Froese et al. (2017), the default prior for relative biomass, B/k was set to 0.2-0.6 (medium) for the start year (B start /k). From FishBase (Froese and Pauly, 2019), the resilience of L. litulon is "Low, " corresponding to a prior r range of 0.05-0.5 (Froese et al., 2017). The prior range for an unexploited population size or carrying capacity (k) was calculated as: where, k low and k high are the lower and upper boundary priors for k for low and high levels of biomass at the end of the time series, respectively, max(C) is the maximum catch value of the time series of catch data, and r low and r high are the lower and upper boundary priors of the r value.
Applying the values of r and k, and the relative biomass for the first year of the time-series of catch data (B start ), the Monte Carlo method was used to filter out suitable r-k pairs. Basic biomass dynamics are described by the following formula: where, B t+1 is the exploited biomass in the subsequent year t+1, B t is the current biomass, and C t is the catch in year t.
(3) is replaced by the following: where, the term 4 B t k assumes a linear decline in recruitment below half of the biomass capable of producing MSY. More detailed equations and concepts can be found in Froese et al. (2017). Given that the assessment was based on just 8 years of catch and effort data, sensitivity analyses were also conducted to investigate the potential effect of the relative biomass prior (B start /k; 0.2-0.6) on estimates of B/B MSY and F/F MSY , compared with a higher prior (0.4-0.8) and a lower prior (0.1-0.4).
The time series of catch and CPUE are displayed in Supplementary Table 4, and the R code for the BSM method, as well as Supplementary Materials describing the method in detail, can be downloaded from R Core Team (2018).

Abundance Maximum Sustainable Yield Method
The AMSY model was used to determine the maximum sustainable yield from the abundance data. The advantage of this method is that it uses only CPUE or other relative abundance time series to evaluate the exploitation pattern and stock status without catch data (Froese et al., 2020). The required input data for the AMSY method included a time-series of CPUE and prior ranges for r and relative stock size B t k in a given year. The timeseries of CPUE data are given in Supplementary Table 4, and the method of obtaining the prior range of r was the same as described above (Froese et al., 2020). Priors for B k together with population dynamics and the proportional factor of upper and lower limits were used to put the observed CPUE into a preliminary MSY framework, which was then refined by Monte Carlo filtering. The AMSY method estimates the relative catch based on several transformations of the Schaefer model, which requires biomass data for two consecutive years (Froese et al., 2020); hence, the relative catch (Catch/MSY), fishing pressure (F), and exploitation level (F/F MSY ) were estimated up to the secondto-last year in the time series. Froese et al. (2020) provide a detailed description of the theory and equations behind AMSY. In our contribution, sensitivity analyses similar to the BSM method were also conducted for AMSY.

Resource Distribution
The spatial distribution of L. litulon abundance in Shandong coastal waters from 2012 to 2019 (except for 2017) is shown in Figure 2. Lophius litulon was caught in 70 fishing areas. The distribution features of this stock in the coastal waters of Shandong in recent years were as follows: (1) the species was widely distributed in the Yellow Sea, but was less common in the Bohai Sea; (2)    *** Indicates extremely significant differences (P < 0.001); ** indicates a significant correlation at the level of 0.01 (P < 0.01); * indicates a significant correlation at the level of 0.05 (P < 0.05). Figure 3 shows the monthly distribution of L. litulon, which was caught in 43 fishing areas. The species appeared more concentrated offshore from January to September and moved west toward the coast from October to December. The highest CPUE was usually in the central and southern coastal areas of Shandong, peaking in November and December. In January, April, May, and October, the CPUE values were low over the entire study area, with the lowest value occurring in October.

Effects of Environmental Factors
A GAM was used to fit the relationship between the relative abundance of L. litulon and environmental factors. To obtain the optimal model, impact factors were screened according to the minimum Akaike Information Criterion (AIC). The full screening process is presented in Supplementary Table 6. The final expression of the model is as follows: Deviation analysis indicated that the cumulative deviation of the selected modeling factors was 51.70% (R = 0.446, P < 0.05). In addition to the influence of month (time factor), the F-tests showed that the environmental factor with the greatest influence was depth (D), followed by latitude (Lat), BT, BS, and longitude (Lon , Table 1). There was an extremely significant correlation between BT and the relative abundance of L. litulon (P = 0.000728). Figure 4 shows the effect of environmental factors, as simulated by the GAM, on the abundance of L. litulon.   (-0.197-4.26) N; (3) L. litulon was densely distributed in the waters with BT varying between 5.8 • C and 10.6 • C; (4) as BS increased, L. litulon abundance decreased, especially in the range from 31 to 33.2 ; (5) there was a positive association between L. litulon abundance and Lon.

Resource Evaluation
The results and confidence intervals for r, k, and MSY evaluated using the BSM and AMSY methods were similar ( Table 2). The exploitation rate (F/F MSY ) was consistently estimated to be > 1.0, and the relative biomass (B/B MSY ) in the last year was predicted to be < 1.0, which implied an unhealthy L. litulon stock status in Shandong coastal waters. Figure 5A shows 19,027 feasible r-k pairs screened by BSM (black dots; gray dots for CMSY), and the darker gray cross was the best r-k pair with a 95% confidence interval (CI) (the light gray cross for CMSY). Figure 5B depicts 5003 feasible r-k pairs identified by AMSY (black dots), and the cross represents the best pair with its 95% CI. Figure 5C (estimated by the BSM method), indicates that L. litulon catches fluctuated greatly over time, with the lowest value occurring in 2013 and the highest in 2018. Figure 5D shows the time-series of CPUE data overlaid with the estimated biomass that would achieve MSY (gray line). The CPUE predicted by AMSY also exhibited significant fluctuations between 2012 and 2019, with the lowest value in 2017 and the highest in 2018.
The relative biomass trajectory (B/B MSY ) produced by BSM is depicted in Figure 6A, and the gray area represents its 97.5% CI. B/B MSY has been declining since 2015, and the B/B MSY value in the last year of the time series was 0.538 (< 1). Figure 6B shows that the AMSY model estimated that the relative biomass trajectory (B/B MSY ) decreased gradually before 2014, stabilized between 2014 and 2016, and increased gradually after 2016, but the end point of B/B MSY was 0.801 (< 1), which was closer to the level MSY. Figure 6C shows that the exploitation rate (F/F MSY ) curve predicted by BSM, increased gradually after 2013, and the terminal F/F MSY value was 2.41 (> 1). Figure 6D shows that the exploitation rate estimated by AMSY gradually decreased after 2014, and the F/F MSY in the second to last year was 1.47 (> 1).
The results of the sensitivity analyses are shown in Figure 7 and Supplementary Table 7, which indicate that the change in priors had a high impact on the outcome of the BSM, being very low in the AMSY model. Herein, selecting the appropriate prior values is crucial to avoid bias in the BSM estimates.

Spatial-Temporal Dynamics and Environmental Influences
In this study, the marine area with the highest frequency of L. litulon occurrence from 2012 to 2019 was the central Yellow Sea (34.5 • -37 • N, 121 • -124 • E).  investigated the relative resource density and distribution of L. litulon from 1985 to 2009 and found that the area with the highest density of L. litulon in the Yellow Sea was between 34 • -35 • N and 122 • -123 • E. In this study, although the mean CPUE in some years (e.g., 2014 and 2019) was lower than that in adjacent years, L. litulon was more widely distributed than in the 1985-2009 period.
During the 9 months of this study, L. litulon was widely distributed over the entire marine area from February to April, but the CPUE was low. This was possibly because spring is the spawning season for L. litulon, and either the catches comprise mostly recruits with light individual weights or the individuals are too small to be caught by nets. Conversely, the high CPUE in autumn (from September to November) might be attributed to biomass accumulated during the summer fishing moratorium. Li et al. (2012) found that ST had a significant effect on L. litulon CPUE in the southern Yellow Sea in spring.  also pointed out that the catch yield of L. litulon in the central and southern parts of the Yellow Sea was significantly correlated with the ST. This study showed that BT was the main environmental factor significantly affecting the abundance of L. litulon in the Yellow Sea. Furthermore, L. litulon was concentrated in areas of the Yellow Sea and Bohai Sea with low temperatures, high salinities, and depths >18 m. Wang et al. (2013) showed that L. litulon in Haizhou Bay and its adjacent sea areas (within the scope of this study) had a strong negative correlation with temperature, salinity, and depth in winter, but was positively related to surface pH. Herein, we suggest that further studies should consider more environmental factors to improve the accuracy of the correlation estimates. In addition, different resource-modeling approaches should be developed and adopted.

Stock Status and Suggestions for Sustainable Development
In many developing countries, including China, most of the species caught on a large scale have not yet been assessed, nor do they have an adaptive management plan to guarantee their sustainable use and protection (Costello et al., 2012). Stock assessments are generally lacking for fish that gradually become a dominant part of the catch and progressively increase their economic value, as is the case for L. litulon. In 2002, the total annual catch of L. litulon in the East China Sea was ranked ninth of 153 species, accounting for 1.14% of the total catch, and its existing stock size was estimated to exceed 2,000 tons (Lin and Zheng, 2004). The frequency and relative resource density of L. litulon in the central and southern Yellow Sea showed an increase from 1985 to 2009 . However, owing to the strong fishing pressure in the south of the Yellow Sea, the caught individuals of L. litulon became smaller and the population structure became younger .
In this study, the results of the BSM and AMSY approaches showed that, although there were some differences in precision and biomass trajectories between different methods, the same conclusions were evident. That is, the population of L. litulon is now in an unhealthy and overfished state, and the biomass remains below the level that can produce MSY. Our results are consistent with the findings of Wang et al. (2020), indicating that the resource status of L. litulon is unacceptable in the coastal waters off Shandong.
Interestingly, the BSM showed that the status of L. litulon stock is worsening, while the AMSY model suggested it is improving. Although AMSY uses the surplus yield model of the filtered r-k pairs to predict catches that conform to the CPUE time series and the priors, its estimates of the exploitation index (F/F MSY ) are normally given with wide margins of uncertainty (Froese et al., 2020). The accuracy and applicability of AMSY will be affected by how closely the stock abundance and catch follow the assumptions of the surplus yield model, so it may be less suitable than BSM for management purposes. Nevertheless, AMSY should be well suited for estimating the productivity index (r) and relative stock size (B/B MSY ), so it may be useful in the management of data-poor stocks. Furthermore, sensitivity analysis demonstrated that the BSM was more sensitive to the initial relative biomass prior than AMSY. In other words, B/B MSY and F/F MSY estimates generated by the BSM model were greatly affected by the selection of priors, which might be attributed to the short time series of catch and CPUE used in this study. In such cases, reasonable prior ranges are important for obtaining reliable estimates.
Since the early 1990s, China has adopted a series of fisheries management systems to develop and protect its marine resources, including fishing moratoria, fishing permits, fishing effort controls, quota systems, and resource allocation. However, declines in marine fishery resources continue to occur. Therefore, cautious management of L. litulon resources in coastal waters of China is required.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because our data are sourced from existing fisheries data.

AUTHOR CONTRIBUTIONS
ZZ completed the data collection and analysis and wrote the first draft. YW provided guidance on research methods and writing ideas and participated in data processing. YW, CL, and WX conceived and designed the study, and revised the first draft. SL provided the data. All authors contributed to the manuscript and approved the submitted version.