Changes in Distribution Patterns for Larimichthys polyactis in Response to Multiple Pressures in the Bohai Sea Over the Past Four Decades

Understanding patterns of change in the distribution of species among their critical habitats is important for analyzing population dynamics and adaptive responses to environmental shifts. We investigated spatio-temporal changes in small yellow croaker (Larimichthys polyactis) using eight alternative models fitted to data from bottom trawl surveys conducted in the Bohai Sea each spring (spawning period) and summer during 1982–2018. These models included different combinations of local sea temperature, fishing pressure, and individual climate index (i.e., North Pacific index, NPI, and West Pacific index, WPI) as explanatory variables. Selection of the most parsimonious model for each season was based on Akaike’s Information Criterion (AIC). The model with NPI as its only explanatory variable was used as a base case for pre-analysis. In spring, a spatio-temporal model with sea temperature as a quadratic effect, plus the spatially varying effects of a climate index and fishing pressure was selected, as the AIC value of this model was reduced by 41.491 compared to the base case model without these effects. In the summer after spawning, the spatio-temporal model with WPI as a climate index covariate lagged by 1-year best explained the spatio-temporal distribution patterns of the stock. The results suggested that small yellow croaker populations significantly decreased in biomass in the Bohai Sea over the study period. A statistically significant northeastward shift in the center of gravity (COG) and a contraction in the distribution range occurred in summer throughout the study period (p<0.05). During the spring sequence (1993–2018), a statistically significant northeastward shift in the COG was also found (p<0.05). Our results showed that biomass-density hotspots of small yellow croaker in both seasons have shrunk or disappeared in recent years. Overall, these findings suggest that the spatio-temporal patterns of the populations in their spawning, feeding and nursery grounds have been influenced over the past 40 years by multiple pressures, and population density in the southwestern areas of the Bohai Sea declined faster and more drastically than in the northeastern areas. This study has important implications for developing targeted spatial conservation measures for small yellow croaker at various stages of its life history under different levels of stress.

Understanding patterns of change in the distribution of species among their critical habitats is important for analyzing population dynamics and adaptive responses to environmental shifts. We investigated spatio-temporal changes in small yellow croaker (Larimichthys polyactis) using eight alternative models fitted to data from bottom trawl surveys conducted in the Bohai Sea each spring (spawning period) and summer during 1982-2018. These models included different combinations of local sea temperature, fishing pressure, and individual climate index (i.e., North Pacific index, NPI, and West Pacific index, WPI) as explanatory variables. Selection of the most parsimonious model for each season was based on Akaike's Information Criterion (AIC). The model with NPI as its only explanatory variable was used as a base case for pre-analysis. In spring, a spatiotemporal model with sea temperature as a quadratic effect, plus the spatially varying effects of a climate index and fishing pressure was selected, as the AIC value of this model was reduced by 41.491 compared to the base case model without these effects. In the summer after spawning, the spatio-temporal model with WPI as a climate index covariate lagged by 1-year best explained the spatio-temporal distribution patterns of the stock. The results suggested that small yellow croaker populations significantly decreased in biomass in the Bohai Sea over the study period. A statistically significant northeastward shift in the center of gravity (COG) and a contraction in the distribution range occurred in summer throughout the study period (p<0.05). During the spring sequence , a statistically significant northeastward shift in the COG was also found (p<0.05). Our results showed that biomass-density hotspots of small yellow croaker in both seasons have shrunk or disappeared in recent years. Overall, these findings suggest that the spatio-temporal patterns of the populations in their spawning, feeding and nursery grounds have been influenced over the past 40 years by multiple pressures, and population density in the southwestern areas of the Bohai Sea declined faster and

INTRODUCTION
Key ontogenetic habitats for fish populations (spawning grounds, nursery grounds, etc.) support their reproduction and continual replenishment, and are often areas subject to intensive human activities and high levels of exploitation. In the context of global environment change (increasingly intensified impacts of human activities exacerbating adverse climatic conditions), understanding the response mechanisms to environmental change among fish habitats is prerequisite for sustainable development of marine fisheries and protection of their supporting ecosystems (Lotze et al., 2006;Browman and Skiftesvik, 2014;Jin et al., 2015). Studies of sustainable fish production as marine ecosystems evolve under environmental stressors is at the forefront of contemporary fisheries research .
Tracking and understanding distribution shifts of fish stocks and patterns of expansion or contraction in their spawning, nursery, and feeding grounds is fundamental to effective fisheries management. This understanding can help fisheries scientists and managers to clarify the response mechanisms of stocks to environmental changes (Astarloa et al., 2021) and the downstream effects of these changes on the delivery of fisheries ecosystem services such as food production (Myers and Worm, 2003;Lotze et al., 2006;Cheung et al., 2013). In so doing, it facilitates the establishment of a scientific theoretical framework for ecosystem-based adaptive, and in some instances restorative, management of fisheries (Botsford et al., 1997;Pikitch et al., 2004;Beddington et al., 2007;Jin et al., 2015;Thorson, 2019a).
There is increasing evidence that anthropogenic pressures (e.g., fishing) have led to dramatic distribution changes in the critical habitats of many marine species populations (e.g., Liu et al., 1990;Engelhard et al., 2014;Han et al., 2021;Grüss and Thorson, 2019). Fishing can greatly reduce stock biomass and alter population structure (Li et al., 2012;Bell et al., 2015). It can also give rise to changes in life history characteristics and phenotypic evolution (Sun et al., 2022), typically leading to shrinkage or displacement of the area occupied by exploited fish populations (Engelhard et al., 2014;Bell et al., 2015). Seawater pollution caused by human activities can also seriously damage the functions of key habitats. Species diversity and abundance of fish eggs in polluted waters can become significantly reduced, leading to a narrowing and spatial shift in spawning grounds (Cui et al., 2003).
Environmental stressors (e.g., temperature changes, climate changes) can also lead to spatial shifts among fish stocks (e.g., Blanchard et al., 2005;Pinsky et al., 2013). Many studies have already shown that temperature change is the key factor that initiates migration between key habitats (Feng and Yang, 1955;Liu et al., 1990;Jin et al., 2005), and is also an important environmental determinant of stock distribution. Disruption of migration and restriction of habitable environments may make it difficult to obtain adequate catches from existing fishing grounds in the future (e.g., Overholtz et al., 2011;Cheung et al., 2013;Bell et al., 2015;Su et al., 2015). Water temperature also affects the survival, development and hatching of fish eggs (Bian et al., 2014), thereby causing fluctuations in population replenishment and changes in fish distribution. In some cases, species distributional shifts are related to the nature of climate change rather than increasing temperature per se (Pinsky et al., 2013). A study in the Bay of Biscay showed that the regional climate index provided a better explanation than local environmental variables (SST and chlorophyll) in the spatio-temporal distribution patterns of dolphins (Astarloa et al., 2021), although local effects may be more influential on the trophodynamics of teleosts than mammals.
Environmental stressors can also have profound indirect effects on the early life stages of marine fish in both nursery and feeding grounds. It has been found that long-term changes in the water temperature of the North Sea cod (Gadus morhua) habitat caused by global climate change has altered species composition, population structure and abundance of copepods that the cod larvae rely upon for food. High larval mortalities from their inadequate food supply, have consequently impaired stock recruitment giving rise to dramatic fluctuations and an overall decline in North Sea cod stocks over the last few decades (Beaugrand et al., 2003;Richardson et al., 2004). In many instances, depending on their specific life history stages, carnivorous fish populations are far less directly sensitive to temperature and other stressors than their planktivorous prey. Ignoring indirect effects from declines in prey may render exploited fish species without adequate protection and preclude their optimal utilization.
Critical habitats are increasingly recognized by researchers as being affected by combinations of multiple stressors. For example, climate warming and fishing pressure have combined to influence a century of distributional shifts of Gadus morhua in the North Sea (Engelhard et al., 2014). In China, multiple stressors have caused shifts in ecosystem structure and function during the spawning season in Laizhou Bay (Jin et al., 2013). Jin et al. (2013) showed that top-down effects have been the dominant influence on species composition over the last half century due to increased fishing pressure, whilst bottom-up effects have increased over the last thirty years due to strong changes in environmental stressors (other stressors such as temperature, nutrients, and salinity). These findings imply that fisheries scientists and managers need to identify and monitor the multiple stressors affecting fish populations to effectively manage exploited fish stocks.
The Bohai Sea is a semi-enclosed shallow sea in China's warm temperate zone and is connected with the Yellow Sea to its east (Figure 1), which is a key habitat integrating spawning, nursery and fishing grounds (Liu et al., 1990). Many highly migratory fishes of the Yellow and Bohai Seas with high socio-economic importance spawn and forage in the Bohai Sea, with the latter playing a crucial role in replenishing most fish populations in the Yellow Sea Large Marine ecosystem (Jin et al., 2005;Jin et al., 2015;Bian et al., 2018). Small yellow croaker (Larimichthys polyactis, northern Yellow Sea-Bohai Sea stock) is the most iconic and important among exploited species in this region, relying on habitats in the Bohai Sea for its spawning and nursery phases of life before overwintering in the fishing grounds of the Yellow Sea (Liu et al., 1990). Despite the important implications of the relationship between small yellow croaker and multiple stressors during its spawning and nursery periods for stock conservation, scant research has been conducted on this topic. Therefore, there is an urgent need to better elucidate the patterns of variation in the spatial distribution of the stock in its spawning, nursery and feeding grounds of the Bohai Sea, how these patterns change over time, and how these patterns respond to variations in anthropogenic and environmental stressors. This will help with developing strategies for sustainably managing the stock. In the present study we formulated a spatio-temporal variable coefficient model and applied it to the Bohai Sea for the first time to identify the combination of stressors affecting longterm distributional changes in small yellow croaker.

MATERIAL AND METHODS
In this study, fishery independent catch rates (kg.km -²) of small yellow croaker from the Bohai Sea (northern Yellow Sea-Bohai Sea stock) during spring and summer from 1982-2018 were fitted using a spatio-temporal delta-Gamma model. The main objective was to identify shifts in distribution patterns indicative of range expansion or contraction of the stock in its spawning, nursery, and feeding grounds under multiple environmental stressors. Initially, we developed eight alternative spatiotemporal models of the stock for each season, constructed from combinations of explanatory variables that included the effects of sea temperature, climate index, and fishing pressure. The Akaike Information Criterion (AIC; Akaike, 1974) was used to select the most parsimonious from among the eight alternative models in each instance. We included a set of regional climate indices, representing the main modalities of changing conditions in the North Pacific Ocean shown to have important effects on fish populations in the Bohai and Yellow Seas (e. g. Liu et al., 2017). Using a pre-analysis of each season, the best climate index was selected as a model covariate, described below in more detail. Next, we used estimates from the eight models to identify the relative importance of local temperature, climate index, and fishing pressure for each season. Finally, changes in centers of gravity (COGs) i.e., centroid of stock abundance or biomass, and effective area occupied by the stock were analyzed using the model with the smallest AIC in each season, to capture the distributional shifts in geographic range of small yellow croaker in its spawning, nursery and feeding grounds and its relationship with multiple stressors. In this study, the spatially varying coefficient model (SVC, Thorson, 2019b) was used to represent the effects of annual indices (i.e., climate index and fishing pressure index) in our spatio-temporal model.

Data Used in This Study
The data for small yellow croaker in the Bohai Sea analyzed in the present study were all from scientific i.e. fishery independent, surveys using bottom trawl at fixed stations (Figure 1), conducted by the Yellow Sea Fisheries Research Institute (YSFRI). The spring (May) survey years for which records were available included 1982, 1993, 1998, 2004, 2011Summer (August) survey years included 1982, 1992-1998, 2009-2010, 2012, 2015. Local commercial fishing vessels were hired as the research vessels for trawling (1-hour tows at a speed of 3 knots). Except for 2011 (mono trawl), the surveys were conducted using tandem commercial fishing vessels (pair-trawling). The power of the survey vessels before 2011, 2011 and after 2011 was 200 horsepower, 350 horsepower and 490 horsepower, respectively. Detailed information about these fishing vessels can be found in some published literature (Zhang et al., 2015;Shan and Jin, 2016;Li et al., 2018). The net mouth had a width of 22.6m, height of 6m, circumference of 109.62m, and cod-end mesh size of 20mm. Samples from all stations were identified to the level of species or lowest taxon as far as possible. Biomass was estimated as weight, abundance as number of individuals from each shot.
The environmental stressor data in this study included local sea temperature, climate index, and fishing pressure index for the period 1981-2018. Remotely-sensed sea temperature data were used because temperature data were not collected for the complete survey time series and the entire Bohai Sea area. Daily sea surface temperature (SST) data for the study area from May 1981-2018 were downloaded from the Optimum Sea Surface Temperature (OISST) Database (https://www.ncei.noaa. gov/products/optimum-interpolation-sst) of the National Centers for Environmental Information, NATIONAL OCEANIC AND ATMOSPHERIC ADMINISTRATION (NOAA). The spatial accuracy of the data is 0.25°× 0.25°, which matches the precision of our research grid. We converted the Daily SST into an estimate of mean SST (May) for each survey. Considering that the water depth of the Bohai Sea is relatively shallow, with an average of only 18.7m, and the seawater temperature is low with vertical mixing during spring producing a relatively uniform temperature throughout the water column, SST was used as a proxy for the temperature of each layer of the water column, including the bottom temperature (Jin et al., 2005;Radlinski et al., 2013). The station-measured bottom and surface temperature data obtained during the spring survey in 2017 showed a strong correlation (r = 0.98). Therefore, it is reasonable to assume that SST in spring may have a potential relationship with population dynamics of small yellow croaker in the relatively shallow waters of the Bohai Sea even though it is a benthopelagic species. In summer, SST cannot be used as a substitute for bottom temperature because the upper and lower layers of the water column are not sufficiently mixed, and the 2017 survey data only showed a weak correlation between SST and bottom temperature in August (r = 0.6). Therefore, we downloaded the bottom temperature data with coarse spatial accuracy (0.5°× 0.5°) for the study area from August 1981-2018 and interpolated them. This bottom temperature data information was available at the National Marine Data Center of China (http://mds.nmdis.org. cn/pages/dataViewDetail.html?dataSetId=81).
Pacific Decadal Oscillation (PDO, Mantua et al., 1997), North Pacific index (NPI, Trenberth and Hurrell,1994), West Pacific index (WPI, Barnston and Livezey, 1987), North Atlantic Oscillation (NAO, Hurrell, 1995;Hurrell and Deser, 2010;Hurrell et al., 2003), Northern Oscillation Index (NOI, Schwing et al., 2002), Niño 3.4 index (Trenberth, 1997), Southern Oscillation Index (SOI, Trenberth, 1984) and Arctic Oscillation index (AOI, Thompson and Wallace, 1998;Thompson and Wallace, 2001) data (downloaded from https:// psl.noaa.gov/data/climateindices/list/) and their 1-year advance time series that considered the delayed/lag effect (represented by -_lag1, such as PDO_lag1.) were utilized as climatic pressure indices in the present study. These indices were downloaded for spring months (March-May), summer months (June-August), and the annual individual climatic index values for each study season (spring and summer) used in the present study were the average values of the individual indices for the seasonal month groups for each year. Oceanographic indices summarize the physical status of marine ecosystems (Grimmer, 1963;Kidson, 1975), established by many studies to have direct or indirect effects on marine fishery populations (Mantua and Hare, 2002;O'Leary et al., 2018). These climate indicators are briefly described and summarized in the appendix (Supplementary Table 1). In this study, time series of these regional indices and their respective advance time series that incorporated a 1year lag effect were chosen as climate pressure indices because they are well recorded and updated in the Bohai Sea, and many studies have shown them to be largely associated with changes in the North Pacific marine ecosystem (Nakata and Hidaka, 2003;Overland et al., 2008;Tian et al., 2008;Shen, 2012;Tian et al., 2014;Liu et al., 2017;Yuan et al., 2017;He, 2002;Liu et al., 2021). Some of these indices have been shown to be strongly correlated with yield changes in small yellow croaker .
Finally, fishing mortality has been recognized as a main driver behind variation in distribution in the Bohai and Yellow Seas (Xu et al., 2003;Lin et al., 2016). No accurate commercial catch and effort data were available for the small yellow croaker stock in the Bohai Sea. Fishing power estimates were derived from records provided by the Bureau of Fisheries and Fishery Administration of Ministry of Agriculture and Rural Affairs (Ministry of Agriculture), China. The fishing power for a given year during the period 1982-2018 was obtained as the logarithm of the total engine power of the commercil fishing vessels from the Shandong, Liaoning, and Hebei provinces of China and one Chinese city (Tianjin). In this study, the total vessel power (excluding the vessels operating in the Yellow Sea) was calibrated using internal information obtained by the YSFRI from the fisheries management departments of these provinces. However, with the development of technology and the improving social economy, the number and proportion of high-powered vessels operating in the fishery are increasing, which means that the growth rate in fishing pressure is faster than that of the total fishing vessel's power. Fishing vessel power records had been divided into three categories for the period 2003-2018 (total power above 441kW, between 45 and 440kW, and below 44kW); and into five ranges prior to 2003, but the thresholds of five ranges before 1993 differed from those during 1993-2002. Therefore, we converted the proportion of highpowered fishing vessels into a calibration coefficient for fishing capacity enhancement. The implementation of a summer fishing moratorium reduces the time available for fishing, which affects changes in fishing pressure, so we also needed to account for variations in the proportion of fishing time expended from year to year. Therefore, we used the product of fishing power × calibration coefficient for fishing capacity enhancement × proportion of fishing months for the period 1982-2018 as the fishing pressure index (FI): Where Power All, y is the total power of fishing vessels in year y; Power MG, y is the total fishing vessel power in the maximum fishing vessel power group in y; Maxrange LL,y is the lower limit of the maximum fishing vessel power grouping threshold for year y; Months a, y is the number of months available for fishing in y. We also compared the performance of the two fishing pressure indicators namely the logarithm of the index FI and the logarithm of the total power of fishing vessels Power All, y in the model (see the model description below) developed in this study. Our analysis showed that the logarithm of the index FI was better than the logarithm of the Power All, y , so we used the former as a proxy of fishing pressure (Supplementary Table 2).

Spatio-Temporal Modeling
We developed eight alternative spatio-temporal models for the small yellow croaker stock of the Bohai Sea for the spring and summer seasons. These models variously included or excluded the effects of local temperature, fishing pressure and climate pressure indices (e.g., NPI), with the most parsimonious model chosen as the optimal one for each season based on the lowest AIC (Akaike, 1974). These models were a set of vector autoregressive spatio-temporal models, more precisely, delta-Gamma generalized linear mixed models (delta-Gamma GLMMs) with or without the option of adding a spatially varying coefficient (SVC) model, to account for spatiotemporal structures (unmeasured changes that vary over time) and spatial structure (unmeasured changes in long-term patterns, i.e., stable over time) at a fine scale Thorson, 2019a;Thorson, 2019b). The models also accounted for differences in relative fishing efficiency introduced by the different survey vessels . We fitted these spatio-temporal models to the catch rate data of small yellow croaker collected during the survey in each season (spring and summer). Specifically, the following eight models comprising up to 3 covariates were fitted, M1: a model with no covariates; M2: a model with the quadratic effect of local temperature, representing a dome-shaped response to local temperatures; M3: a model with the spatially-varying effects of the individual climate indices (e.g., NPI) represented using an SVC model option (Thorson, 2019b), and which was used for pre-analysis, with the climate index in the model with the lowest AIC being selected for each season; M4: a model with the spatially-varying effect of fishing pressure; M5: a model with the quadratic effect of local temperature and the spatially-varying effect of fishing pressure; M6: a model with the quadratic effect of local temperature and the spatially-varying effect of climate index (this climate index was selected from the M3 pre-analysis); M7: a model with the spatially-varying effect of fishing pressure and the spatially-varying effect of climate index (r this climate index is selected by M3 pre-analysis); and M8: a model with the quadratic effect of local temperature, the spatially-varying effect of fishing pressure, and the spatially-varying effect of climate index (this climate index was selected from the M3 pre-analysis). Because some of the survey catch rate data was zero-inflated, we chose the spatio-temporal delta-Gamma GLMM  which comprised a binomial-GLMM and a Gamma-GLMM. This GLMM combination fitted the 0/1 (encounter/ non-encounter) data and the data of the stations where small yellow croaker were caught (non-zero catch data), respectively, and multiplied the predicted values of the two parts to obtain the final biomass-density estimates (Lo et al., 1992;Grüss et al., 2019b;Han et al., 2021). In this study, the spatio-temporal delta-GLMM was employed not only to analyse spatio-temporal distribution patterns of small yellow croaker density in spring and summer, but also to identify the extent to which COGs and the effective area occupied by the stock may have shifted during 1982-2018, as described below.
Since model M8 contains the largest number of covariables (including SST, climate index and fishing pressure), its structure is described here; the structures of the M1-M7 models are similar. The binomial-GLMM with a logit link function and linear predictors was used to estimate the probability of small yellow croaker being caught (encounter) p i at sampling station s (i).The binomial-GLMM included two Gaussian Markov random fields which represent respectively the spatial variation and the spatio-temporal variation in encounter probability: is the relative fishing efficiency for the v th survey vessel at site s(i) in year t(i); w (p) s(i) is the spatially correlated variability in encounter probability at the station s(i) where sample i; ϵ s(i),t(i),1 is the spatially-varying effect of the climate indices (e.g., NPI) on encounter probability at s(i) in t(i); and x (p) s(i),t(i), 2 is the spatially-varying effect of fishing pressure on encounter probability at station s(i) in year t(i). Prior to being used in the models, both the T and T² covariates were standardized to have a mean of zero and a variance of one; this transformation implied that g (p) 1 T (p) and g (p) 2 T 2(p) had a standard deviation equal to g (p) 1 and g (p) a multivariate normal distribution and, in the case of the ϵ (p) t , temporal variation was expected to conform to a random-walk process (Thorson et al., 2016a) in time to ensure that year intervals or/and stations that have not been sampled do not present mean-reversion: Where s 2 pr is the variance of vessel effects (encounters);R (k) is the correlation among stations as a function of decorrelation distance k;s 2 pw and s 2 pϵ are the pointwise variance of w (p) theand ϵ (p) t , respectively; P t,1 is the climate indices (e.g., NPI) and q s,1 P t,1 is the climate indices (e.g., NPI) effect; and s 2 px, 1 is the pointwise variance of the climate indices (e.g., NPI) effect.P t,2 is the fishing pressure; and q s,2 P t,2 is the effect of fishing pressure and s 2 px, 2 is its pointwise variance. The calculation method of the R terms can be found in Thorson et al. (2015).
Similarly, the non-zero catch rate of small yellow croaker r i at station s i was estimated by the Gamma components (has a log link function) and linear predictors, including two Gaussian Markov random fields accounting for spatial variation and spatio-temporal variation in non-zero catch rate, respectively: where most of the parameters of Eq. (4) have the same meaning as the parameters of Eq.
(2), except that they apply to log-catch rate.
To improve computational efficiency, following Han et al. (2021), 100 "knots" n j = 100 were specified to approximate all the w (p) , w (r) , ϵ (p) t ,and ϵ (r) t , over a fixed spatial domain W, so that the value of each variation term is tracked at each knot (Shelton et al., 2014;Thorson et al., 2015). These knots were distributed over a grid (Figure 2) developed for this study. All t were tracked at each knot by the model, and the value of each term at a prediction grid cell is interpolated from the value of three knots surrounding that grid cell (see Grüss et al., 2020a;Grüss et al., 2020b;Grüss et al., 2020c). The 100 knots provided a compromise between the estimation accuracy and run time of the GLMM model. As the number of knots increased, the parameter estimates and predictions became similar.
According to the results obtained by multiplying the predicted values of the binomial-GLMM and the predicted values of the Gamma-GLMM, we mapped the biomass-density of small yellow croaker in spring and summer separately on the Bohai Sea prediction grid. Next, the biomass of the small yellow croaker stock in year t,B t was estimated, as: j,t,2 are random effects .
To measure the shifts in spatial distribution of the stock over time, we calculated its annual COGs in spring and summer during 1982-2018. The eastward COG, X t , was calculated by (Thorson et al., 2016a;Thorson et al., 2016b;Thorson and Barnett, 2017): where x j is the easting value (in km) in knot j. The northward COG, Y t , can be shown in a similar way, except that x j is replaced with y j ,the northing value (in km) in knot j, in Eq. (6).
To determine the annual range expansion/contraction of the stock in spring and summer, we estimated the annual effective area occupied by the stock during 1982-2018. Effective area occupied is shown by the ratio of predicted biomass (Eq. (5)) over average density, D t , which is calculated by (Thorson et al., 2016a;Thorson et al., 2016b): The models were implemented using the "VAST" package (https://github.com/James-Thorson-NOAA/VAST Thorson, 2019a) within the R statistical software environment on Windows (R Core Team, 2021). The estimation of fixed effects, random effects and their standard deviations, as well as the standard deviations of derived quantities were mainly calculated by the Laplace approximation, the stochastic partial differential equation method (Lindgren et al., 2011), and the generalized delta method (Kass and Steffey, 1989) implemented using the R package "TMB" , and the bias-correction estimator developed in Thorson and Kristensen (2016). Their detailed computational descriptions are documented in Thorson et al. (2015) and Thorson (2019a). The convergence test of the spatio-temporal model comprised two parts: (1) the gradient of marginal logarithmic likelihood of all fixed effects should be less than 0.0001; and (2) the Hessian matrix of the second derivative of negative logarithmic likelihood must be positive definite. In addition, the "VAST" package also has a set of diagnostics, which can further confirm whether the spatio-temporal models were reasonable, i.e., that (1) the observed encounter frequencies were within the 95% prediction interval of the predicted encounter probability; (2) the diagnostics for the positive-catch-rate component using a q-q plot showed the residuals lay along the one-to-one line; and (3) mapping Pearson residuals for encounter-probability and positive catch rates of the small yellow croaker stock showed that there were no visible patterns in the Bohai Sea.

Understanding the Relative Importance of Three Stressor Covariables in Explaining Distribution Patterns of Small Yellow Croaker
We employed the eight models (with various combinations of the effects of local temperature, climate indices, and fishing pressure) for each season, and applied the method used in Thorson (2015) to analyze the relative importance of multiple stressors explaining the probability patterns of encounter and non-zero density (estimated by the delta-Gamma GLMM) of small yellow croaker in the Bohai Sea in spring and summer, respectively. Based on the method used in Thorson (2015), we compared the variances of the spatio-temporal variations (the sum of the spatial variation term and of the spatio-temporal variation terms) in the probability of encounters and non-zero density for the eight alternative delta-Gamma GLMMs for each season in this study to determine whether some reduction of the variances caused by local temperature and/or climate index and/ or fishing pressure covariable(s) was/were included in a model. The desired goal of the inclusion of covariables in the models was to minimize the spatio-temporal variability, which represents the latent (unmeasured) variation in the encounter probability or non-zero density probability of the stock Han et al., 2021).  We also performed trend analyses of the local temperature, climate index and fishing pressure index time series in spring and summer as an aid for interpreting the estimates of the most parsimonious (i.e., the smallest AIC value) spatiotemporal model. These trend analyses were performed by the regime shift detection method based on a sequential t-test to detect possible regime shifts, which requires the data to be processed with "pre-whitening" prior to use (Rodionov, 2004;Rodionov and Overland, 2005;Rodionov, 2006).

RESULTS
The Relative Importance of the Fishing Pressure, SST and the Climate Index in Explaining Distribution Patterns In this study, eight models (delta-Gamma GLMMs) for the small yellow croaker stock in the Bohai Sea during spring and summer, were developed (M1-M8) by incorporating different combinations of the effects of fishing pressure, local temperature, and climate indices (e.g., NPI). All the models satisfied tests of convergence and diagnostic requirements (The diagnostic results of the AIC-selected models are shown in Supplementary Figure 1, 2).
Comparisons among the M1-M8 models for each season allowed us to identify the relative importance of fishing pressure, local temperature and the climate indices (e.g., NPI) in explaining the probability patterns of encounter and density of the stock during each specific season. The results from the pre-analysis of M3 in spring indicated that M3 with the NPI covariable had a smaller AIC value than M3 with the other climate index covariables. In spring, we found that when interpreting the probability pattern of encounter, the variances of M1-M8 were all equal to or extremely close to 0, and NPI and fishing pressure were relatively more important than local temperature ( Table 1). Including NPI or fishing pressure in the model, the variance of spatio-temporal variation of encounter probability remained unchanged relative to M1, while including local temperature in the model caused a negligible increase of the variance. Conversely, we found that SST was more important than NPI and fishing pressure in interpreting the non-zero density pattern during spring ( Table 1). The inclusion of local temperature in the model caused a large decrease in the spatio-temporal variance in density, whereas the inclusion of NPI or fishing pressure caused a moderate increase in the variance.
The results of pre-analysis of M3 in summer indicated that M3 with the WPI_lag1 (WPI lagged by 1-year) covariable had a smaller AIC value than M3 with other climate index covariables. We found that WPI_lag1 was more important than local The variation of encounter probability means the sum of the spatial variation terms and spatio-temporal variation terms calculated by the binomial component of the spatio-temporal delta-Gamma generalized linear mixed model (GLMM). The variation in non-zero density probability means the sum of the spatial variation terms and spatio-temporal variation terms calculated by the Gamma component of the model. Percent change in variance, Percent change in variance compared to M1 (no covariates included); SST, sea surface temperature; BT, sea bottom temperature; NPI, North Pacific index; WPI, West Pacific index; WPI_lag1, West Pacific Index time series with 1-year advance (Consider the 1-year delay/lag effect); AOI, Arctic Oscillation index; Fishing pressure, fishing pressure index.
temperature and fishing pressure when interpreting the probability pattern of encounter in summer (Table 1).
Including WPI_lag1 in the model, the spatio-temporal variation in encounter probability showed a large decrease relative to M1, while including local temperature or fishing pressure in the model caused a negligible or small decrease in the variance. Moreover, the WPI_lag1was more important than local temperature and fishing pressure in explaining the probability patterns of density; and local temperature was more important than fishing pressure ( Table 1). Including the WPI_lag1 or local temperature in the model caused a large decrease in spatio-temporal variation in density probability, whereas including fishing pressure in the model caused a large increase in the variance.

Change Patterns of Spatiotemporal Distribution and Range of Small Yellow Croaker in the Bohai Sea
In spring, the AIC value of the spatial-temporal model with the quadratic effect of local temperature, the spatially-varying effect of NPI and the spatially-varying effect of fishing pressure was the lowest ( Table 2) t term in equation (4)) reflected the biomass change, which also declined markedly in 1982-1998; increased largely in 1998-2011; decreased sharply in 2011-2014; and increased modestly in 2014-2018. Since 2011, biomass density in spring decreased significantly, accompanied by the disappearance of the hotspot area (Figure 3 and Supplementary Figure 3). This disappearance/contraction was more pronounced in the western (Bohai Bay) and southern (Laizhou Bay) areas of the Bohai Sea than in the mid-eastern regions. After 2011, the density of the stock in the Yellow River Estuary and the southwest areas of Laizhou Bay was close to zero and showed a slow recovery; especially in the southwest areas of Laizhou Bay, which had a negligible recovery in 2017-2018. In 2018, the central area of Liaodong Bay close to the Liaodong Peninsula remained the only density hotspot in the Bohai Sea (Supplementary Figure 3).
In summer, the spatio-temporal model (M3) with the spatiallyvarying effect of WPI_lag1 had the lowest AIC ( Table 2). This AICselected model M3 predicted that the hotspots of small yellow croaker in the Bohai Sea would be found in the area near the line connecting the junction of Bohai Bay and Laizhou Bay to the northern area of Liaodong Bay in the summer of 1982-2018 ( Figure 4). We also found that the predicted hotspots moved left (west) or right (east) along the diagonal between years. M3 also predicted that the biomass decreased markedly in general in the summer between 1982 and 2018, leading to substantial distribution changes of the stock in their nursery and feeding grounds over the entire study period (Figure 4 and Supplementary Figure 4). The b (r) t terms also showed a significant overall decline during 1982-2018. Since 2009, the density of the stock in the Bohai Sea in summer decreased significantly, accompanied by shrinkage of the high-density areas (Figure 4 and Supplementary Figure 4). The shrinkage was more obvious in Bohai Bay, the eastern area of Laizhou Bay and the central area of the Bohai Sea than in Liaodong Bay. In 2012 and 2016, the density in the Yellow River estuary remained relatively high in the area south of 39°N. By the second  Figure 4).
Combining the trends of the stock biomass density in the Bohai Sea in the two seasons, we found that the biomass: (1) increased extremely slowly in the early 1980s to early 1990s; (2) subsequently decreased significantly during the early 1990s to late 1990s; (3) increased again largely in the early 2000s to late 2000s; (4) then decreased sharply in the early 2010s to mid-2010s; and finally (5) exhibited an extremely slow recovery trend after 2017. In spring, the COGs predicted by the AIC-selected M8 suggested that the COG of the stock had moved southward and westward in spring of 1982-1993, and then northward and eastward in spring between 1993 and 2018 ( Figure 5). Over the period 1982-2018, changes in both the COGs in spring were non-significant (both p > 0.05, a two-sided Wald test was used for all significance tests). However, the spring variations in eastward and northward COGs during 1993-2018 were significant (p=0.015 and 0.034, respectively). In summer, the COGs estimated by the AIC-selected M3 showed that the COG variation trend was relatively stable in the summer of 1982-2009; and then the COG moved largely to the north and east over the period 2009-2018 ( Figure 5). Between 1982 and 2018, eastward changes in the COG in summer (i.e., the large move of the COG of the stock towards the east of Bohai Sea) were statistically significant (p = 0.0004, a two-sided Wald test was used for all significance tests). In addition, the northward COG changes of small yellow croaker in summer (i.e., the large move of the COG of the stock towards the north of Bohai Sea) were also significant (p = 0.0024) from 1982 to 2018.
In spring, the effective area occupied of the stock predicted by the M8 indicated that range contraction had occurred during the study period from 1993 to 2018 ( Figure 5). However, from 1982 to 2018, the effective area occupied by the stock in spring did not change significantly (p >0.05). In summer, the effective area occupied in the Bohai Sea predicted by the M3 showed that the range of the stock was reduced ( Figure 5). Statistically, the change in effective area occupied by the stock in summer from 1982 to 2018 was also significant (p = 0.046).

Trend Analysis of Local Temperature, Climate Index and Fishing Pressure Index Time Series
From 1982 to 2018, the local temperature anomalies in May showed an oscillatory rising trend, and the mean SST difference between the former and the latter five years was 1.96℃. The SST time series showed that regime shifts in temperature occurred in 1995/1996 and 2013/2014 respectively ( Figure 6A), mirroring the switch from a southward to a northward shift of the small yellow croaker COG in spring in 1993 and in 2014 respectively. The time series indicated ( Figure 6B) that there were no regime shifts in the spring NPI during the period 1981-2018. The change in trend of the cumulative sum of NPI anomalies in Only projections for years where sampling for the stock were available (i.e., 1982, 1993, 1998, 2004. 2011, 2014 and 2016-2018) are presented.  , 1982, 1992, 1998, 2009-2010, 2012-2013 and 2015-2018) are presented. spring was, to a certain extent, similar to the trend in the COG in spring.
In summer, there was a large seasonal oscillation in local temperature anomalies during August of each year ( Figure 7A).

DISCUSSION
To understand its ecological response to multiple stressors during the spawning and nursery phases of its life history, we developed a spatio-temporal model for the small yellow croaker's populations in the Bohai Sea in spring and summer during 1982-2018. Ours is the first spatio-temporal implementation of the SVC model to determine the effects of annual indices (climate index and fishing pressure index) in the Bohai Sea. Many studies have established that spatio-temporal models which account for both spatial and spatio-temporal variation can produce more accurate estimates, thereby providing increased reliability in the scientific basis for decision-making about how best to manage marine ecosystems, their dependent fish populations, and the fisheries those populations sustain (Rassweiler et al., 2014;Thorson et al., 2015;Grüss et al., 2019b;Han et al., 2019;Thorson, 2019a;Thorson, 2019b;Han et al., 2021).
In the spatio-temporal model we developed, spatial variation refers to variation in biomass density that for the small yellow croaker in the Bohai Sea is assumed to remain constant over time. Spatio-temporal variation refers to the interannual variation in biomass density of the stock. Spatial variation and spatio-temporal variation represent, respectively, the fundamental ecological niche of the stock in the Bohai Sea and the response of the stock to the unmeasured environment stress (Grüss et al., 2017;Thorson, 2019a). Spatio-temporal models can also be used to obtain information on the stressors that cause ecological responses in fish populations and to analyze the relative magnitude of the effects of these stressors (i.e., analyzing the relative importance of the covariates) by including covariates to improve the proportion of response variation explained by a model (Thorson, 2015;Grüss et al., 2020a;Han et al., 2021). Based on these characteristics of the spatio-temporal model, this study identified multiple pressure factors influencing the population dynamics of small yellow croaker during spawning, nursery and feeding phases of their life history in the Bohai Sea.
The model results in spring and summer were expected to accord with those from a study in the East Bering Sea (Thorson, 2019b), which showed that the use of SVC models, including temperature or/and annual index (cold pool), produced more parsimonious models that better described the stock being modeled. The relationship between fishing pressure and local sea temperature with the distribution of small yellow croaker in summer was not as great as that in spring. We suspect that implementation of the summer fishing moratorium system meant that the stock dynamics during summer was not directly affected by fishing. Han et al. (2021) argued that the inclusion of local covariables and/or annual index in a similar model would not be appropriate to describe the stock being modeled in many regions; and the study of Ducharme-Barth et al. (2022) also adds a supporting case for this view. The present study supports and further develops the idea that such regions are not static, and that the inclusion of local or/and distant regional covariates in spatio-temporal models varies greatly in its suitability for describing the modeled stock in different seasons.
Insights about the spawning grounds and reproductive environment of the Bohai Sea stock of small yellow croaker in this study were predominantly consistent with those from previous studies conducted during the 1950s-1980s (Feng and Yang, 1955;Liu et al., 1990), with some additional points. Previous studies indicated that small yellow croaker entered the Bohai Sea from the Yellow Sea and divided into two routes FIGURE 5 | Eastward center of and northward center of gravity (COG; in km) and effective area occupied for the stock in spring and summer of 1982 −2018, predicted by the AIC-selected model. The 95% confidence intervals are represented as shaded areas. Only projections for the years where sampling for the stock were available (i.e., 1982, 1993, 1998, 2004. 2011, 2014 and 2016-2018 for spring; 1982, 1992, 1998, 2009-2010, 2012-2013 and 2015-2018 for summer) are given.
in late April: the north route entered Liaodong Bay in mid-May; the south route continued to branch, one branch reached the spawning grounds of the Yellow Estuary of Laizhou Bay in early May, and the other branch reached the northern area of Bohai Bay (38°30'N-39°00'N, 118°00'E-118°30'E) in mid-May. Our results showed that even in mid to late May, a significant percentage of the stock remained distributed in the central Bohai Sea. Our model results also indicated that the range of spawning grounds in Bohai Bay was larger than determined in previous studies, and mainly concentrated in the northwestern area of Bohai Bay; and that in Laizhou Bay, in addition to the area around the Yellow River estuary, the stock was also distributed in the southern area of Laizhou Bay during the spawning period. The spawning grounds are crucial for replenishing the populations in the entire Yellow Sea Large Marine Ecosystem, because they attract many spawning adults (Jin et al., 2005). Our model results imply that there is a need to take more measures to focus on protecting the spawning grounds in Bohai and Laizhou Bays. Therefore, this study provides fundamental support for marine spatial planning to ensure spawning grounds are sufficiently protected to maintain sustainable fisheries.
Our model indicated that the main feeding grounds (predicted hotspots) were in the area near the line connecting the junction of Bohai Bay and Laizhou Bay to the northern area of Liaodong Bay in summer (August). This location was consistent with the main distribution area of the August fishery yield of the stock in the Bohai Sea observed during 1971-1972(Xu and Chen, 2009, Figure 2 of the paper). However, the common feeding ground in Laizhou and Bohai Bays differed from what was identified in a previous study by Lin (1991), which determined that the common feeding ground of the stocks in Laizhou and Bohai Bays was located in the nearshore waters of Bohai Bay in northwestern Shandong Province. In summer, M3 indicated a statistically significant shift to the northeast and range contraction of the stock in their feeding and nursery grounds over the period 2009-2018. This finding mirrored a decrease in the recruitment capacity of small yellow croaker in the Bohai Sea after 2009, with its northeastern area declining relatively slower than the other regions of the Bohai Sea. As per the spring results, this finding also implied a decline in ecological function (cradle) in Bohai and Laizhou Bays in recent years and the need to focus on restoring their functionality.
Next, we discuss the role of local temperature, climate index and fishing pressure in driving the decadal spatio-temporal patterns of the stock in its spawning, nursery and feeding grounds in the Bohai Sea.
Previous studies have established that fishing pressure is an important driving force for the distribution patterns of marine fish stocks (Jørgensen et al., 2008;Bell et al., 2015). In the Yellow and Bohai Seas, fishing pressure and its heterogeneous distribution have also been identified as important factors leading to changes in fish distribution (Xu et al., 2003;Li, 2011;Wang et al., 2012;Lin et al., 2016). In a spatio-temporal model, the total power of fishing vessels is far from the best method to represent fishing pressure (Han et al., 2021), because it ignores the exponential increase in fishing efficiency of individual vessels and the changes in available fishing periods caused by management regulations. In this study, we recognized this issue, and in the fishing pressure index we constructed, we accounted for the reduction in fishing periods caused by the summer fishing moratorium and for the rapid increase in fishing efficiency arising from the increase in the proportion of high-powered fishing vessels (which tend to have more sophisticated technical equipment). Although there are many other factors affecting fishing pressure, we can only narrow the distance between this index and the actual fishing pressure based on the limited information available, but we believe that the index we used is more effective and realistic than the total power of fishing vessels, and is the most effective proxy of fishing pressure in this area at present. We therefore fitted a spatio-temporal model that included this fishing pressure index effect, modeling it as an annual index. In spring, based on the AIC comparison of the eight alternative spatio-temporal delta-Gamma models, we found that the effect of fishing pressure was second only to the effect of SST. Fishing pressure effect was also included in the final chosen model, M8. We believe that it is possible that heterogeneity in the distribution of fishing pressure has caused a more pronounced decrease in the density of small yellow croaker in the western (Bohai Bay) and southern (Laizhou Bay) areas of the Bohai Sea than elsewhere. This is evident as a decrease in the effective area occupied and a northeasterly shift in the COG of the stock as predicted by the spatio-temporal model. In summer, fishing pressure did not have a direct impact on the small yellow croaker in August for most of the inter-annual time series due to the closure to commercial fishing during summer since 1995 (after several adjustments later, i.e., the duration of the moratorium was increased). The spatio-temporal distribution of the stock in August was indirectly affected by fishing pressure via its influence on the spawning parental stock in spring. This explains why in this study, we found that the model including the fishing pressure effect reduced AIC compared with the model without this covariate, but that its AIC was not the smallest.
Previous studies have shown that fishing pressure has a substantial effect on the biomass of small yellow croaker (Liu et al., 1990;Lin et al., 2016). This is confirmed by the results from applying the spatio-temporal models in spring and summer in our study. In terms of trends in the status of the fish stocks in the Bohai and Yellow Seas, Liu et al. (1990) showed that serious overfishing occurred from the early 1970s to the late 1980s. Our results indicate that a substantial part of the large increase in the biomass from the early 2000s to the end of the 2000s was associated with a steep decrease in fishing pressure during that period (Figure 8). Fishing pressure dropped to its lowest in 2009, providing an important respite for the replenishment of stocks, and was directly reflected in the high biomass observed in the following summer of 2009 (no data were available in the spring of 2009, but survey data for the autumn of 2009 also showed that the density of the stock was much higher than in previous years), and indirectly led to the high biomass density of adult small yellow croaker observed during the spring of 2011. This was partly supported by the changes in the density of small yellow croaker in the overwintering grounds of the Yellow Sea (Han et al., 2021), where the density in the northern area (considered to be the overwintering grounds of the Bohai-Northern Yellow Sea stock) and the central overwintering grounds (considered to be the overwintering grounds of the Bohai-Northern Yellow Sea and the central Yellow Sea stocks) increased in 2008 relative to previous years, and relative to the southeastern overwintering grounds of the Yellow Sea. Unfortunately, although the total power of fishing vessels remained stable after 2009, the proportion of high-powered vessels sharply increased (implying faster speed; increased of fish holding capacity; increased periods of continuous operation at sea; advanced fish finding sonar technology, etc.), leading to a sharp increase of fishing pressure, which in turn led to a deterioration in the status of the small yellow croaker resource. Although the summer season closure was extended for another month from 2017, the resource did not recover significantly. The need to constrain fishing pressure has gained more and more attention from fisheries managers, with the Ministry of Agriculture and Rural Affairs of the People's Republic of China (2017) issuing the "Notice on Further Strengthening the Control of Domestic Fishing Vessels and Implementing the Management of Total Marine Catch", which specified clear control targets and an implementation schedule for vessel and catch limitation. In combination these imposed a strong fishing pressure control measure. Our study also shows that reduced fishing pressure can indeed promote the recovery of small yellow croaker stocks. In addition, our results also highlight the destructive effect of highpowered fishing vessels on stocks. At the same time, measures need to be taken to ensure that there are enough spawning adults of small yellow croaker during May and enough numbers of juveniles leaving the Bohai Sea in October to begin their overwintering migrations.
Model selection indicated the importance of incorporating local temperature into the models of the stock during the spawning season ( Table 2). SST may explain a large proportion of the spatio-temporal variation of non-zero density probability ( Table 1). This finding is consistent with previous results that changes in ocean temperature affect the onset of spawning migrations and the distribution of spawning grounds (bottom temperature 10−13°C) for the stock (Liu et al., 1990).
During the summer feeding and nursery periods, we found that local temperature also explained a large proportion of the spatio-temporal variation of non-zero density probability and a smaller proportion of the spatio-temporal variation of encounter probability in August, but the proportion explained was less than the climate index WPI_lag1. Model selection also indicated that the spatio-temporal model with a local temperature effect was less suitable for the modeled systems than the spatio-temporal model with WPI_lag1 effects. This is strongly related to the wide temperature range (generally 14-26°C) to which small yellow croaker are subjected during their feeding period (Zhao et al., 1987). WPI_lag1, on the other hand, contains other information in addition to temperature, which will be discussed later. As found in Astarloa et al. (2021), the summer local environmental variable SST in our study was unable or rarely able to capture the association between the environmental and ecological processes because of the time lag of species' response and the inherent nonlinear nature of stock dynamics (Hallett et al., 2004). We hypothesized that the distribution of the stock during the summer feeding period may be strongly correlated with the distribution of its primary prey, and the time-lagged driving factors for the primary prey distribution over a period of time in advance can be used as time-lagged covariates among the explanatory effects in the summer spatio-temporal models.
Many species require non-local and regional mechanisms to explain their spatio-temporal distribution shifts and density changes (HilleRisLambers et al., 2013;Heino et al., 2017). Even highly compressed annual indices (climate/oceanographic indices) that lose spatial differences in distribution sometimes contain additional information about species distributions compared to local temperatures (Thorson, 2019b). Thorson (2019b) concluded that the SVC (spatially varying coefficient) for an annual oceanographic/climate index represents a flexible way to proxy the habitat selection as fish respond differently among years in accordance with their changing physical environments. This is the case for the summer model results in this study, where additional information on the distribution of small yellow croaker during the summer feeding period, compared to local temperature, was included in the regional covariate, namely annual climate index WPI_lag1. WPI_lag1 represented the combined effect of lagged and geographically distant habitat conditions, which influence the choice of feeding grounds. WPI_lag1 was the best predictor (lowest AIC) for explaining the biomass density of small yellow croaker in summer. The results showed that biomass density and encounter of small yellow croaker in summer were negatively correlated with WP_lag1 (and, positively correlated with WPI). The east-west and south-north movement of the East Asian jet stream leads to positive and negative WPI patterns, so WPI is related to various aspects of the East Asian climate. WPI mainly affects precipitation, temperature and occurrence of tropical cyclones in East Asia (Wallace and Gutzler, 1981;Barnston and Livezey, 1987;Choi and Moon, 2012). We also found a symmetric relationship between WPI and local temperature during 1982-2018, with positive WPI anomalies often corresponding to negative SST anomalies (Figure 7). Astarloa et al. (2021) suggested that climate indices often affect the abundance of high trophic level predators by regulating their food resources rather than directly affecting them. Therefore, we could assume a potential bottom-up process, in which WPI_lag1 influences the small yellow croaker by regulating the effects of precipitation and other environmental effects on prey.
The results of the spring model showed that the climate index NPI also contains a small part of this combined effect, but does not explain the distribution of small yellow croaker to the extent that the climate index WPI_lag1 does in summer. The model results showed that NPI in spring was negatively correlated with the encounter and density of the stock. Spring NPI is thought to affect the temperature of the Yellow and Bohai Seas to promote gonad development of small yellow croaker , which could promote the spawning migration of the stock from the overwintering grounds to the spawning grounds of the Bohai and Yellow Seas. In addition, the spawning grounds are also affected by runoff into the sea, and spring NPI may affect the runoff from the River (Jiang et al., 2008;Zhi et al., 2015).
In conclusion, local sea temperature, climate index and fishing pressure exerted measurable effects on the spatial distribution pattern (fundamental ecological niche), shifts in that pattern, and range expansion/reduction of small yellow croaker populations during spawning and feeding seasons.
There are still many aspects of these patterns that cannot be explained by these three indices ( Table 1, Figures 5). The fractions of these patterns that were unexplained by the three indices, were to a large extent driven by the decreases in small yellow croaker biomass in 1982-2018 (Figures 3, 4), and partially by other unconsidered factors (shoreline reclamation, seawater pollution) or unknown sources. Jin et al. (2015) summarized the impediments to the sustainable output of Chinese inshore fishery resources, especially in the Yellow and Bohai Seas, and concluded that ecological disasters such as red tides and the fragmentation or functional loss of spawning and nursery grounds caused by high-intensity human activities such as large-scale reclamation projects for coastal development, landbased pollution and mariculture have seriously impaired the replenishment and sustainability of fish stocks. Therefore, it is necessary to develop meaningful proxy indices that can appropriately represent the impacts of these high-intensity human activities for future research and fisheries assessments.
This study has contributed to our further understanding of the spatial distribution pattern of small yellow croaker and its relationship with multiple pressures in spawning, post-spawning feeding and nursery grounds in the Bohai Sea during summer. Other studies of the spatial distribution pattern of small yellow croaker habitat in the Bohai Sea were based on data from the 1950s to 1980s (Feng and Yang, 1955;Liu et al., 1990;Lin, 1991), and therefore do not provide the necessary insight into the present distribution pattern to enable responses to the critical needs of spatial management and protection in the Bohai Sea, as the cradle for the fishery in the Yellow Sea Large Marine Ecosystem. The results of this study are generally compatible with those of previous studies, but our research has highlighted that the patterns of spatial distribution of the small yellow croaker stock in their spawning and feeding grounds in the Bohai Sea have changed significantly over the past 40 years. In addition to the changes in spatial pattern identified in this paper, adaptive responses among the reproductive characteristics of the stock to environmental changes have also been identified by Zeng et al. (2005) who showed that individual fertility of small yellow croaker with the same body length had increased significantly in 2004 compared with 1964. Since the 1980s, the age of sexual maturity of small yellow croaker has significantly advanced, and the reproductive stock is mainly composed of recently matured fish (Li, 2011). We therefore suggest that spatio-temporal models should be further developed to accommodate the effects of these changes in reproductive characteristics and explore their relative importance in explaining patterns in spatial distribution shifts and range expansion/contraction. In this way, future resource management of the Bohai Sea and even the whole Yellow Sea Large Marine Ecosystem can include development of special conservation plans for the restoration of fish stocks based on previous and anticipated environmental and fishing patterns (Grüss et al., 2018;Grüss et al., 2019a;Grüss et al., 2014).
This study has provided important information for the management of small yellow croaker resources in spawning, nursery and feeding grounds in Bohai Sea. It has also provided a spatio-temporal model framework for the study of other highly migratory species in the Bohai Sea, thereby supporting spatial conservation plans and ecosystem-based fisheries management measures. Importantly, this study has demonstrated that SVCs are suitable as an annual index which serves as an important planning tool by including not only climatological information (Thorson, 2019b) but also fishing pressure into species distribution models (SDMs). It also illustrates that other highintensity human activity pressures (such as pollution and reclamation) could be incorporated into SDMs in the future. This will promote the application of indices of fishing pressure and other highly impactful human activities for predicting changes in the spatio-temporal distributions of fish stocks via SDMs, and ultimately promote the implementation of ecosystem-based fishery management.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Ethics Committee, Yellow Sea Fisheries Research Institute.

AUTHOR CONTRIBUTIONS
QH and XS conceived this study. QH conducted to methodology, analysis, writingoriginal draft, visualization. XJ made substantial contribution to data curation, writingreview & editing. HG made substantial contribution to conceptualization, writingreview & editing. CS made contribution to writingreview & editing. All authors contributed to the article and approved the submitted version.