Impact of Climate Change on Wintering Ground of Japanese Anchovy (Engraulis japonicus) Using Marine Geospatial Statistics

The distribution and fluctuations in abundance of small pelagic species such as anchovy are largely affected by climate change. We hypothesized that the future projected rise in temperature will result to a northward shift of Japanese anchovy (Engraulis japonicus) habitat and a subsequent increase in relative abundance. To test this hypothesis, we explored the link between Japanese anchovy abundance and environmental conditions using machine-learning and statistical models. The models were fitted with catch per unit effort (CPUE) as the response variable and remotely sensed data of sea surface temperature (SST), sea surface chlorophyll-a (Chl-a), assimilated information of sea surface salinity (SSS), meridional and zonal ocean currents, and depth as environmental covariates. Our results showed that the abundance of E. japonicus was significantly influenced by environmental factors. In particular, salinity front and SST highlight strong relationships with winter CPUE distribution. Based on these models, the results reinforced our hypothesis and showed that the warming ocean will drive a substantial shift in Japanese anchovy habitat in the China seas. SST and CPUE showed negative correlations with the El Niño Southern Oscillation (ENSO) index. These findings underpin ramifications of the climate-driven habitat shift of small pelagic fish species on the regional marine ecosystem in the China seas.


INTRODUCTION
Climate change is considered as one of the major drivers of the fluctuations in abundance and distributions of small pelagic fish resources such as anchovy (Lecomte et al., 2004). For instance, European anchovy in the northern British waters is expected to shift northward in response to climate change (Pinnegar et al., 2017). In the China seas, Japanese anchovy is also highly susceptible to climate change (Ma et al., 2019a). One of the climate phenomena that impacts fish resources is the El Niño Southern Oscillation (ENSO), which is a coupled atmosphere-ocean variability in the tropical Pacific Ocean. This event is also defined by indices such as the Niño1 + 2, which better captures the ENSO signal in winter (Trenberth and Stepaniak, 2001;Hanley et al., 2003;Gen et al., 2010;Yan et al., 2017). During an El Niño, the surface temperature changes in the Yellow Sea (Wei et al., 2010), and the increase in precipitation enhances the riverine discharge from the Yangtze River, thus leading to salinity changes (Park et al., 2015). Many studies have found that fishery resource variations are closely linked to ENSO events. For instance, ENSO is shown to affect the hotspots and distributions of chub mackerel in the East China Sea (Yasuda et al., 2014) and also causes regional changes in the density of young Bluefin tuna in the East China Sea (Kitagawa et al., 2006). Fluctuation in abundance between the sardine and anchovy is influenced by ENSO (Brierley and Kingsford, 2009). Few studies to date have explored the relationship between climate change and the distribution of Japanese anchovy, specifically in the climate-sensitive waters of the China seas where ENSO was shown to affect its early life history (Tian et al., 2004;Kim et al., 2005). Thus, the impact of present and future marine environmental changes to fish resources deserves our attention. We proposed a hypothesis that a climate-driven temperature rise will result in the northward shift of anchovy habitat and increase of anchovy relative abundance. We tested this hypothesis using the climate projections from the IPCC (Intergovernmental Panel on Climate Change)-RCP (Representative Concentration Pathway) scenarios. The IPCC forecasts future climate change and discusses the impacts of climate change, future risks, and adaptation strategies. The R of the IPCC project global mean surface temperatures describe four scenarios (RCP 2.6, RCP4.5, RCP6.0, and RCP8.5) based on different greenhouse gas emissions (Nurdin et al., 2017;Silva et al., 2019). By the end of the 21st century (2081-2100), the IPCC-RCP scenarios projected that the global surface temperature would rise between 1.0 • C and 3.7 • C compared to the 1986-2005 baseline (IPCC, 2014b). Thus, this raises the urgent need to study the impacts of future climate change on fisheries. Using the Japanese anchovy as a model species to study this topic is likely to disclose processes and mechanisms similar to other small pelagic fisheries across the globe.
Japanese anchovy (Engraulis japonicus) is a widely distributed species and is particularly abundant in the northwestern Pacific region. It is a small pelagic fish that undergoes diurnal vertical migration and is considered a key species linking lower and upper trophic components of the food web (Zhao et al., 2008). Japanese anchovy is a forage fish that moves between its spawning, feeding, and wintering grounds (Figure 1). The separate stocks of Japanese anchovy do not exist in the Yellow Sea and East China Sea (Zheng et al., 2015). During the onset of its spawning migration in late March, Japanese anchovy migrates north and northwest of the basin and spawns from June to August. After spawning, E. japonicus commences its feeding migration and shifts south from November toward its wintering grounds in the Yellow Sea and the East China Sea, where it has the most concentrated and stable stock in winter (December to March) (Zhao, 2006). Therefore, our analyses are focused on the wintering ground of Japanese anchovy in these waters, and it exists in the eastern and southern regions of the Yellow Sea. It is also concentrated in the northern East China Sea and around FIGURE 1 | Map showing the Japanese anchovy migration route modified from Zhu (1990) and Zhao (2006). The different regions comprising the China seas are indicated on the map, and overlain are the migration route (purple arrows), spawning (red), feeding (blue), wintering (black), and winter fishing (green) grounds.
Jeju Island, with the most considerable wintering ground in the central waters of the southern Yellow Sea (Kwon et al., 2012;Huang et al., 2014;Niu et al., 2014). We delineated the regions of the eastern and southern Yellow Sea and northern East China Sea as the winter fishing grounds (WFG). The monthly averaged CPUE within the winter fishing ground was calculated as the WFG CPUE index (ton h −1 ). Chinese industrial trawl vessels catch Japanese anchovy, and the government implemented the summer fishing moratorium for protecting fishery resources (Jiang et al., 2009). Since the 1990s, the stock biomass of Japanese anchovy in the Yellow Sea has significantly declined from 4.12 million tons in 1993 to 1.75 million tons in 2000 based on acoustic surveys (Zhao et al., 2003). The fishing yield of Japanese anchovy declined from 1.26 million tons in 2001 to 0.52 million tons in 2009, which then increased to 0.95 million tons in 2015; the Japanese anchovy catch accounts for 10.56% concerning the total fish catches in China (Yuan and Zhao, 2016).
In order to effectively monitor and manage the fishery for making sustainable utilization of marine resources, a number of geospatial statistical methods have been used to explore the relationships between environmental factors and species distributions. Among these statistical tools are machine learning methods such as random forest (RF) (Smoliñski and Radtke, 2017), generalized linear models (GLM) , and generalized additive models (GAM) (Murase et al., 2009). As the applicability of models is dependent on the target species and regional differences, we constructed all three models to identify the best suited model for projecting the abundance and distribution of wintering Japanese anchovy in the northwestern Pacific region. Remote sensing technology could provide environmental information of suitable temporal and spatial resolution (Klemas, 2012). Moreover, Geographic Information Systems (GIS) technology provides analytical tools to extract environmental ranges from the spatial surveys and map the essential species habitats (Valavanis et al., 2008). Likewise, Simple Ocean Data Assimilation provides numerical model-derived data for comprehensive and accurate oceanographic information continuously corrected through the assimilation of observational data (Carton and Giese, 2008). It can be used complementary to remote sensing data for species-environment studies. Thus, the combination of GIS, remote sensing technology, and statistical models can provide scientific information for marine resource management (Mugo et al., 2011) and the design of appropriate policies to protect fishery resources.
The objectives of this study are (1) to elucidate wintering ground distributions of Japanese anchovy using multiple models; (2) to analyze the impacts of various environmental factors on Japanese anchovy, and (3) to investigate the effect of projected climate scenarios on its spatial and temporal distributions.

Study Area and Target Species
The China seas consist of the waters of the Bohai Sea (BS), Yellow Sea (YS), East China Sea (ECS), and the South China Sea, and this study was conducted in regions of the Yellow Sea and East China Sea (26-42 • N and 118-128 • E) (Figure 1). The YS is a semi-enclosed shelf between China and the Korean Peninsula while the East China Sea is located south of the Yellow Sea (Dong et al., 2011) and links Chinese, Korean, and Japanese waters. The Yellow Sea and the East China Sea are affected by the tides and monsoon winds and influenced by the intrusion of the Kuroshio current (Xing et al., 2020). The China seas are also among the most productive fishery areas in the word (Ma et al., 2019b), which represent 66.10% of the total catch of the Japanese anchovy in the northern Pacific (FAO, 2017).

Fishery Data
Daily catch data of Japanese anchovy in China seas based on the double trawl surveys were provided by the East China Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences. These data were digitized and compiled into a monthly database from 2009 to 2015 (Table 1), which comprised the daily fishing locations (latitude and longitude), fishing dates, and catch per unit effort (CPUE), expressed in terms of tons/hour-boat day (ton/fishing boat-hours). We preprocessed the fishery data and averaged the CPUE when fishing activities were operated in the same location for many consecutive days within a month. An abnormally high CPUE or an abnormally low CPUE during a short fishing interval at the same fishing location were considered to be outliers. We then matched the fishery information with the environmental data at each specific date to obtain the database for fishing location and environmental data.

Environmental Variables
Remotely sensed Level 3 sea surface temperature (SST) and sea surface chlorophyll-a (Chl-a) were derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data from 2009 to 2015 1 ( Table 2). However, the numerical modelderived sea surface salinity (SSS), meridional (V) and zonal (U) current velocities from Simple Ocean Data Assimilation (SODA) v3.3.1 were downloaded from the APDRC data repository 2 . From the SSS data, salinity fronts were identified based on Cayula and Cornillon single-image edge detection algorithm (Cayula and Cornillon, 1992). Bathymetric data were also downloaded from the online database 3 . All data were resampled into a 0.01 • spatial resolution.

Future Projections and Climate Index
The projection was based on different IPCC-RCP scenarios with temperature increases at 1.0 • C, 1.8 • C, 2.2 • C, 3.7 • C, for RCPs 2.6, 4.5, 6.0, and 8.5, respectively, at the end of the 21st century. The CPUE distributions of Japanese anchovy were predicted based on the temperature increases for different scenarios at the end of the 21st century. The anomalies were calculated as the difference between predicted CPUE based on different scenarios and the predicted CPUE in 2015. The index of Niño1 + 2 was 1 https://oceancolor.gsfc.nasa.gov/ 2 http://apdrc.soest.hawaii.edu/data/data.php 3 https://www.ngdc.noaa.gov/mgg/global/ The data sources and original spatial and temporal data resolutions are noted.
used to represent the ENSO events (Hanley et al., 2003) and was downloaded from the online data repository 4 . The wintering fishing ground CPUE index was used to analyze the relationship with mean SST in China seas and the Niño1 + 2 index.

Model Construct
We constructed a total of 63 GLMs, 31 GAMs, and 6 RF models using the stepwise add variable approach, and the intermediate results were shown in Supplementary Material. Here, environmental data (SSS, Depth, Chl-a, SST, U, and V) were used as the predictor variables in all models to predict the CPUE distribution of the Japanese anchovy. Marine Geospatial Ecology Tool (MGET) of ArcGIS 10.5, an add-in program for applying advanced analytical methods (Roberts et al., 2010), was used to detect the frontal features. Generalized linear models generalizes the traditional linear models, which is an expansion of the maximum likelihood estimation in probit analysis that combines both systematic and random components in the model (Nelder and Wedderburn, 1972). GAM is a semi-parametric extension of a GLM and assumes that the functions are additive and the components are smooth. GAMs can deal with highly non-linear and nonmonotonic relationships between the response and predictor variables (Guisan et al., 2002). Here, GAM is constructed using the "mgcv" package (Wood, 2006). The collinearity of all predictor variables was examined using the variable inflation factor (VIF) (Supplementary Material 1), where factors without collinearity have VIF values less than 5 (Cornic and Rooker, 2018). As all predictor variables have VIF below the cut-off, they were all used for constructing the models. In the construction of GAMs, the smoothing function used was the thin plate regression spline with "shrinkage" to get poorly performing predictor variables. Due to the low model contribution of V, this factor was disregarded to improve the model's predictive performance (Roberts et al., 2010). The CPUE follows a continuous distribution, logarithmic transformation was used to normalize the asymmetrical distribution, and a factor of 0.1 is commonly used in CPUE standardizations (Supplementary Material 2.1). GAM is modeled using the Gaussian distribution family (Solanki et al., 2017). The GAM Equation (1) 4 https://www.ncdc.noaa.gov/teleconnections/enso/indicators/sst/ and the GLM Equation (2) implemented were presented as follows: Where g is the link function, Y is the response variable (CPUE), X i is the predictor variables (SSC, Depth, Chl-a, SST, U, V), f i is a spline smoothing function for each model predictor X i , a, and β are model constants, and ε is a random error (Wood, 2006). Finally, RF model is a classification method that uses multiple decision trees, which are used as statistical classifiers composed of tree predictors. RF generates decision trees based on predictor variables from the original database, and decision trees provide classification by a majority vote from trees (Breiman, 2001). The RF model is constructed using the "Random Forest" package (Liaw and Wiener, 2002).

Model Selection Method
Pearson's correlation coefficient was used to test the linear relationship between two variables by using a t-test; the tests of statistical significance levels used were 0.05 and 0.01 (Liu et al., 2014(Liu et al., , 2018. In GLM and GAM, the addition of different predictor variables leads to changes in the deviance explained and Akaike's Information Criterion (AIC). Model selection was based on the contribution of the predictor variable, the highest values of deviance explained, and the lowest values of AIC (Johnson and Omland, 2004). For the RF models, the chi-square test and AIC results were not computed. Instead, the model selection was only based on the increases in deviance explained. Finally, the databases were used to construct RF, GLM, and GAM for the period between 2009 and 2014, and the best models of RF, GLM, and GAM were used to predict CPUE using monthly raster composites of environmental predictors in 2015. Model performance was evaluated by comparing the predicted and actual fishing CPUE in 2015.

Comparison of RF, GLM, and GAM
The final selected models for GLM, GAM, and RF (model, predictor variable, AIC, p-value, and deviance explained) are shown in Table 3. In GAM, the best-performing model has five predictor variables (SSS, Depth, Chl-a, SST, U) that were highly significant (p < 0.01) with the lowest AIC (217.53) and the highest deviance explained (36.7%). For GLM, the final selected model with three significant predictors (SSS, Chl-a, V; p < 0.01) had the lowest AIC (303.84) and a high explained variance (9.30%). In RF, the addition of predictor variables resulted in an increased deviance explained; hence, the model with six variables (SSS, Depth, Chl-a, SST, U, and V) was the selected model.
The predicted CPUE ranged from 0 to 2 tons per hour (ton h −1 ) using the selected models for RF, GLM, and GAM. A histogram of residuals and normal quantile-quantile plots indicated that the distribution of residuals of GAM adequately conformed to the application of Gaussian distribution (Supplementary Material 2.2). Figure 2 shows GAM was the best model by comparing the relationships between the predicted and actual CPUEs for each of the selected models. The GAMbased model showed the highest correlation (R 2 = 0.72, p < 0.01), which was greater than those from RF (R 2 = 0.19, p < 0.01) and GLM (R 2 = 0.15, p < 0.05).

Preferred Oceanographic Conditions for Japanese Anchovy
Generalized additive models analysis indicated that five variables such as SSS, Depth, Chl-a, SST, and U significantly influenced CPUE of Japanese anchovy (p < 0.01). The importance of each variable was ranked based on the deviance explained and AIC ( Table 4).
The non-linear relationship between environment variables and CPUE is shown in Figure 3. GAM plots can be interpreted as the effect of each environmental variable on the CPUE. A narrow confidence interval represents high relevance. The positive effect of SSS on the CPUE was observed between 28 and 33. From 18.5 to 28 and 33 to 34.5 (Figure 3A), there were negative effects on the fish CPUE and the density of the data points within the former was lower with a wider confidence interval relative to the latter. A positive effect of depth on CPUE was observed between −50 and −110 m, while negative effects were noted from −30 to −50 m ( Figure 3B). For Chl-a ( Figure 3C), a positive effect on the CPUE is evident from 0.71 to 1.25 mg m −3 , and a negative effect was observed between 1.25 and 3.16 mg m −3 . The plot highlighted a positive effect on CPUE beyond 3.16 mg m −3 , but with a wide confidence interval. Considering the confidence interval and the fitted curve, a positive SST effect on CPUE was observed for ranges between 8-14 • C and 19-22.5 • C ( Figure 3D).   The positive effects of U on CPUE were captured between −0.09 and −0.02 m s −1 , between 0.01 and 0.04 m s −1 (Figure 3E).

The Distribution of Japanese Anchovy Prediction Wintering Ground CPUE
The distribution of monthly CPUE (ton h −1 ) of Japanese anchovy predicted during the wintering period (2009-2015) using the GAM model is shown in Figure 4. Japanese anchovy was distributed in the eastern and southern parts of the Yellow Sea and the northern waters of the East China Sea. Its distribution also extended near the Yangtze River estuary and the Jeju Island.

Correlation Between SST and the WFG CPUE Index
The SST showed that it decreased to 9.2 • C in January 2011 and then increased to 10.5 • C in February 2011 (Figure 6), but generally increased in China seas. The WFG CPUE index also tended to have an increasing trend from 2009 to 2015. In December, the WFG CPUE index was highest and declined as SST decreased. SST and WFG CPUE index showed a significant positive relationship (R 2 = 0.69, p < 0.01) ( Figure 7A).

Impact of Environmental Variables on the CPUE of Japanese Anchovy
Our work proposed multiple marine statistical models to elucidate wintering ground distribution of Japanese anchovy, based on remotely sensed data and assimilated information.
In the marine environment, the statistical model has been extensively used to predict species richness, diversity, and geographical abundance distributions (Smoliñski and Radtke, 2017). By comparing the predicted CPUE and the actual CPUE and the predicted and actual fishing ground locations, a multiplemodels approach could generate the robust wintering ground predictions for Japanese anchovy. The spatial distribution of Japanese anchovy cannot be adequately monitored due to technical constraints and economic reasons, and the marine geospatial statistics model as an accurate and cost-effective technology could predict the inevitable gaps and species distribution in the unsurveyed area (Murase et al., 2009). The statistical models show the distribution of anchovy CPUE and infer the wintering ground of anchovy location. The model disclosed the impact of climate change on small pelagic fish by using Japanese anchovy as a model species; the result held the hypothesis that the rising temperature caused by climate change resulted in the shift of anchovy habitat northwardly and the increase of anchovy relative abundance. We could quantitatively evaluate the predictive ability of the marine statistical model through the comparison of different models for fishery resource protection and management application. Freshwater discharge from the Yangtze River affects the salinity of the East China Sea (Delcroix and Murtugudde, 2002). The wide range of salinity preference of Japanese anchovy ( Figure 3A) suggests that it can take advantage of nutrient-rich estuarine waters for foraging. The Changjiang Diluted Water makes the Yangtze River estuary and the adjacent East China Sea increasingly eutrophicated (Chai et al., 2006) and provides excellent feeding opportunities for the Japanese anchovy during their spawning and feeding migration. The concentration of winter fishing locations close to salinity fronts also supports these findings (Figure 8). The wintering migration of the Japanese anchovy is driven by temperature and salinity. A strong salinity front forms west of the Jeju Island in winter (Lie, 1985), where the wintering grounds for Japanese anchovy in the Yellow Sea coincides with the front near the Jeju Island (Huang et al., 2010). The Japanese anchovy catch in the south sea of the Korean peninsula accounted for 25% of the total catch in the adjacent waters (Kim and Kang, 2000). The salinity front created a good feeding environment for Japanese anchovy due to the accumulation of microalgae and zooplankton in the frontal region (Bertrand et al., 2002).
The salinity front is strongest in winter and is associated with topography (Chen, 2009). Depth is also relatively important for the CPUE distribution of the Japanese anchovy, where shallow waters (depth <50 m) generally exerted a negative effect on its CPUE. In the western margin of the North Pacific, depth indirectly plays an important role in the distribution of Japanese anchovy. Japanese anchovy feeds on zooplankton rather than the phytoplankton (Niu et al., 2014). The earlier study has reported that the catch of Japanese anchovy is related to Chl-a concentration (Ihsan et al., 2018). Chl-a and zooplankton were also shown to be positively correlated with the catch of Japanese anchovy (Kim and Kang, 2000). Micro-zooplankton provides food sources for Japanese anchovy during their spawning and feeding migration, and Chl-a plays an indirect role in the distribution of Japanese anchovy. However, there may be a time-delayed relationship between pelagic fish and primary productivity (Yen and Lu, 2016), and Japanese anchovy feed rarely during their winter migration (Zhu, 1990). This may account for low deviance explained obtained for Chl-a. Temperature plays a vital role in Japanese anchovy distribution (Hayashi et al., 2016). The SST range of positive effects on spawning and feeding migration is 19-22.5 • C, and our results show that SST significantly influenced Japanese anchovy CPUE. The range between 8 • C and 15 • C was suitable for wintering Japanese anchovy, and it is concentrated between 11 and 13 • C (Ma, 1989;Li et al., 2007). SST has a positive effect on the abundance of anchovy between 5 • C and 9.5 • C in winter (Niu et al., 2014). Our results were also consistent with the previous researches. The Japanese anchovy migrated vertically in the water layer (Zhao et al., 2008). The temperature experienced by the fish in different depth layers is the variable, and vertical change of temperature may be more influential to the distribution of Japanese anchovy and may explain the low deviance explained by SST.
The salinity front is related to the ocean current in the China seas (Chen, 2009). The current plays an indirect role in the distribution of the CPUE. In the sub-mesoscale physical process, plankton is concentrated as patches, and these patches were related to the spatial distribution of anchovy (Bertrand et al., 2008). Upwelling, that brings nutrient-rich waters to the sub-surface, also promotes primary productivity and copepod production, thus, facilitating the feeding and growth of Japanese anchovy larvae (Nakata et al., 2000).

Impact of Climate Change on Japanese Anchovy CPUE in Winter
Sea surface temperature was associated with biomass of wintering Japanese anchovy (Niu and Wang, 2017), and its catch positively correlates with SST in December (Kim and Kang, 2000). This is also consistent with our findings of an overall increasing SST trend in the China seas corresponding to an increase in WFG CPUE from 2009 to 2015 (Figures 6, 7A). The wintering ground formed in December, and as the SST decreases, Japanese anchovy in the Yellow Sea migrated southward and partially entered the East China Sea (Zhao, 2006). This migration pattern was also captured from our study and was apparent from the respective CPUE decrease and increase in the Yellow Sea and the East China Sea from December to March (Figure 4).
The suitable SST promoted Japanese anchovy migration to northern regions under warm periods over a long-term time scale (Zhou et al., 2015). Our results likewise supported our hypothesis of potential northward habitat shift and CPUE increase for Japanese anchovy under future projected warming (Figure 5). During actual fishing activities, the fishermen carry out fishing operations in areas with high fish density for several consecutive days to ensure a high catch rate (measured by CPUE). If fishing activities cannot maintain a high level of catch rate, then fishermen tend to move to other areas . Thus, the catch rate (CPUE) could be used as an index of abundance (Maunder and Punt, 2004) and is assumed to provide a better indication of abundance than the catch and fishing effort (Lehodey et al., 1998;Chang et al., 2012). It is also effective to use the catch rate (CPUE) to explore the effects of climate change on fish (Syamsuddin et al., 2016). The maximum catch potential in the central and eastern parts of the Yellow Sea is also increasing during 2051-2060compared to 2001(IPCC, 2014a. This suggests that the catch potential of Japanese anchovy  could increase in response to the rising temperatures. Our results showed that CPUE would increase in response to future climatedriven warming, which is also consistent with our hypothesis. As a small pelagic fish, Japanese anchovy exerts a topdown control on zooplankton and bottom-up control on pelagic predators, and hence occupy a central role in the marine food web (Rice, 1995). It is imperative that we pay attention to the migration of small pelagic as it can induce the redistribution of pelagic predators under climate change. The climate impacts on various marine species are increasingly prominent and require robust forecasts of species distribution for conservation and management of fishery resources. Global climate change effects on small pelagics (e.g., anchovies, sardines) are summarized in Table 5. For instance, European anchovy could better adapt to environmental variability. The IPCCbased model predictions for the European anchovy and Atlantic horse mackerel showed potential increase and northward habitat movement (Lenoir et al., 2011;Petitgas et al., 2013;Raybaud et al., 2017). However, the European anchovy stock was expected to decline in southern Europe due to climatic and eutrophication impacts (Macías et al., 2014;George, 2019). California anchovy was also projected to shift northward and likely become a dominant species (Cheung et al., 2015). These earlier results are also consistent with our findings. Moreover, two-thirds of the 36 reported warming-related species shifts in the North Sea highlighted northward habitat movement (Perry et al., 2005). The horse mackerel and anchovy showed a northward expansion of suitable habitats and an increased occurrence probability in the northern British waters by the 2090s (Pinnegar et al., 2017). The copepod abundance and composition will most likely lead to an increase in anchovy abundance under global warming (Stenevik and Sundby, 2007). In the Southern Hemisphere, the habitat distribution of Peruvian anchovy also showed climatedriven poleward shifts (except summer), with decreases in the potential relative in the north and central-south areas (Silva et al., , 2019. Similarly, anchovy and sardine exhibited a large-scale phenomenon of alternation in abundance, associated with climate-associated temperature changes. In particular, the anchovy and sardine abundance showed respective positive and negative correlations against SST, with a stronger correlation between SST and anchovy abundance (Chavez et al., 2003;Checkley et al., 2009;Nakayama et al., 2018). The CPUE and landings of common sardine in the coastal areas off Chile were also projected to decrease in response to climate change (Silva et al., 2015;Yáñez et al., 2016). The biomass and expected profits of the European sardine will obviously decrease with increasing surface temperatures in Iberian-Atlantic (Garza-Gil et al., 2011). However, it would increase in the Bay of Biscay, which is north of the Iberian-Atlantic (Shchepetkin and McWilliams, 2005). The Japanese and Pacific sardines were also expected to expand northward with warming and the latter's suitable habitat was expected to decrease by 50% in the Gulf of California and Mexican Northwest (Okunishi et al., 2012;Cheung et al., 2015;Petatán-Ramírez et al., 2019).
Our result also showed that anchovy CPUE has a significant positive correlation with SST, albeit the Niño1 + 2 showed a significant negative correlation with SST and WFG CPUE index ( Figures 7B,C). During an El Niño (positive Niño 1 + 2), the SST in the China seas was lower than during La Niña (negative Niño 1 + 2) (Jia and Sun, 2002;Wang et al., 2012). One possible mechanism to explain the impact of ENSO on the Japanese anchovy habitat could be provided by looking into the ENSO-related phenomena: ENSO affects the precipitation and subsequent discharge of the Yangtze River that influences temperature and salinity in the China seas. The variability in marine environments affects the wintering grounds and the life history of the Japanese anchovy. For pelagic fish, ENSO is also known to impact the environmental features of their spawning and feeding grounds (Yu et al., 2019). Here, the extent of its impact is likewise evident in the Japanese anchovy wintering ground, when it has the highest annual concentration.

CONCLUSION
The marine geospatial statistical model using environmental predictors provided important information on the wintering ground of Japanese anchovy. A comparison of the three models showed that GAM is the most suitable prediction model for wintering Japanese anchovy. Its distribution was not only affected by SST but also closely related to the salinity front, feeding opportunity, and ocean currents. From SSS images, the salinity front was detected as the main surface feature associated with its wintering ground, highlighting the importance of SSS to the distribution of anchovy. The impact of ENSO on the Japanese anchovy distribution is captured through its influence on the SST on the wintering fishing ground, where CPUE and Niño1 + 2 indices showed a significant negative correlation. We also proposed a hypothesis and explored the influence of future warming on its winter habitat based on four IPCC climate scenarios. Our research reinforced our hypothesis, where temperature increases resulted in a northward habitat shift and CPUE enhancement by the end of the century. This information would be useful in formulating protection policies and promoting sustainable development of fishery resources in the region.

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

AUTHOR CONTRIBUTIONS
All authors listed have made a direct and substantial contribution to this work. SL performed the data analyses and wrote the manuscript. YL conceived the idea for the study and analyses and edited the manuscript. IA edited and proofread the earlier versions of the manuscript. YT collected the data and edited the manuscript. ZY, HY, and JL edited the manuscript. JC provided the fishery data.