Our research focuses on the Food and Agriculture Organization (FAO, https://www.fao.org) Fishing Area 61 within the Northwest Pacific Ocean, spanning from 20°N to 65°N and from 110°E to 175°E. We selected nine highly productive representative species from this region based on FAO publications—Fishery and Aquaculture Statistics (FAO, 2021). Their ecological habitats were determined using information from FishBase (https://www.fishbase.org) and other relevant published sources (Liu and Ning, 2011). These species were distributed in different climatic zones: tropical, subtropical, temperate, and boreal. Furthermore, based on their living depths, they were categorized into pelagic and benthic species. Occurrence data for these species were collected from published literature and two online databases: the Global Biodiversity Information Facility (GBIF, https://www.gbif.org) and the Ocean Biodiversity Information System (OBIS, https://obis.org). To minimize sample bias, we retained one occurrence record per 5.5×5.5 km grid cell (Me et al., 2015). A total of 1552 occurrence records for the nine species were compiled for the modeling process. The ecological characteristics and the number of available occurrence points for each species are detailed in Table 1 . The geographic scope of the study area, along with the distribution points of these species, is illustrated in Figure 1 .

Considering the distribution environments of these species ( Table 1 ) and data availability, we selected ten and twelve predictor variables related to marine surface and benthic environments for pelagic and benthic species, respectively ( Table 2 ). We obtained the environmental data from the online database of Bio-ORACLE v3.0 (https://www.bio-oracle.org). All environmental variables have a spatial resolution of 0.05° (approximately 5.5 × 5.5 km at the equator). These variables were then clipped to the study area using the “crop” function in the R package “raster” for subsequent analyses. Research indicated that ocean temperature, dissolved oxygen concentration, pH, and primary productivity are fundamental controlling factors in marine ecosystems (Fischlin et al., 2007). Sea water velocity, depth, salinity, and slope are commonly utilized as predictor variables in marine biogeography studies (Bradie and Leung, 2017; Assis et al., 2018; Melo-Merino et al., 2020). The mixed layer depth is vital in regulating the exchange of heat and carbon between the atmosphere and the ocean. It supports the majority of primary production within the marine ecosystem, thereby playing a crucial role in sustaining marine life (Sallée et al., 2021). Hence, the mixed layer depth was selected as a predictor variable for pelagic marine species. For benthic species with habitats intrinsically linked to seabed substrates, variables including the topographic position index, terrain ruggedness index, and topographic aspect are deemed significant for their habitats (Pittman and Brown, 2011). To avoid multicollinearity among environmental factors, the two variables with a Pearson correlation coefficient greater than |0.7| are collinear and one of them was removed ( Supplementary Figure 1 ; Dormann et al., 2013; Schickele et al., 2020; Sun et al., 2022). Additionally, we analyzed the collinearity of predictor variables using the Variance Inflation Factor (VIF). The VIF values of the retained variables were all below 3 ( Supplementary Tables 1 , 2 ), indicating the low collinearity among them (Kim, 2019; Sioni et al., 2019; Panzeri et al., 2024). For pelagic species, the dynamic variables of mean ocean temperature (Temp), mean mixed layer depth (Mlay), mean pH (pH), mean salinity (Salt), mean sea water velocity (Cv), mean primary productivity (Phy) in the surface, and the static topographic factors of mean bathymetry (Depth) and topographic slope (Slope) were retained for modeling ( Supplementary Figure 1 ). For benthic species, the dynamic variables of mean ocean temperature (Temp), mean sea water velocity (Cv), mean dissolved oxygen concentration (Do), mean primary productivity (Phy) in the bottom, and the static topographic factors of mean bathymetry (Depth), topographic slope (Slope), topographic position index (Position), and geomorphological features (Aspect) were retained for modeling. To match the temporal scale of the available occurrence data, the environmental data representing current climatic conditions in this study are based on the average environmental conditions of the years from 2000 to 2020.

We categorized the environmental variables into three groups for distinguishing the contributions of different types of variables on species distribution: the physical oceanographic group (POG), the physiology-based group (PBG), and the substrate group (SG) (Pickens et al., 2021). The POG included the predictors of mean bathymetry, mean sea water velocity, and mean mixed layer depth. The PBG included the predictors of mean ocean temperature, mean dissolved oxygen concentration, mean salinity, and mean pH. Mean primary productivity is usually calculated based on phytoplankton abundance, which affects the prey distribution of these important economical species (Roxy et al., 2016; Lewis et al., 2020). Thus, we added this factor to the physiology-based group. The SG included the predictors of topographic slope, topographic position index, and topographic aspect.

To predict and compare the distribution changes of species under future climate scenarios, we accessed future climate data also from the online database Bio-ORACLE v3.0. This database provides free downloads of future environmental conditions under various emission scenarios based on the Shared Socioeconomic Pathways (SSPs) outlined in the Intergovernmental Panel on Climate Change’s (IPCC) Sixth Assessment Report (AR6) (Assis et al., 2024). For our study, we selected two future time periods: 2040-2060 and 2080-2100, and employed three SSP scenarios to represent varying challenges of climate change: SSP1-2.6 (sustainability), SSP2-4.5 (middle of the road), and SSP5-8.5 (fossil-fueled development) (O’Neill et al., 2017). The topographic factors were assumed to be unchanged over time.

We used stacked species distribution models supported by the R package SSDM to model the distribution of nine ecologically significant species in the Northwest Pacific Ocean and their assemblages (Schmitt et al., 2017). Compared to other SDM packages such as ‘sdm’ and ‘biomod2’, the SSDM package enables the concurrent prediction of distributions for multiple species and facilitates assessing their collective assemblage scenarios (Thuiller et al., 2009; Naimi and Araújo, 2016). This package includes nine distinct modeling algorithms: generalized additive models (GAM), generalized linear models (GLM), multivariate adaptive regression splines (MARS), classification tree analysis (CTA), generalized boosted models (GBM), maximum entropy (MAXENT), artificial neural networks (ANN), random forests (RF), and support vector machines (SVM). Due to lack of actual absence data, we used the method based on Barbet-Massin et al. (2012) embedded in the SSDM package to compute and generate pseudo-absence points for each algorithm. This method has been widely used in the SDM researches (e.g., Mainali et al., 2015; Jarnevich et al., 2015; Naimi and Araújo, 2016). Following the method, the regression techniques such as GAM, GLM, and MARS, as well as for MAXENT achieved an optimal model performance when a large number of pseudo-absence points through random selection were generated. However, classification and machine learning techniques, including GBM, CTA, ANN, RF, and SVM, showed a better model performance when a moderate number of pseudo-absence points were generated by constructing geographical exclusion with a 2° buffer. The number of pseudo-absence points for each model per species was showed in the Supplementary Table 3 . The performance of each model was evaluated using a five-fold cross-validation, which was repeated ten times for robustness. This approach randomly partitions the data into a 4:1 ratio, allocating 80% for training and 20% for validation (Guisan et al., 2017). Model transferability and predictive performance were evaluated based on two commonly used evaluation indices: area under the receiver operating characteristic curve (AUC) and true skill statistic (TSS). AUC is calculated based on the retained test data during cross-validation and signals robust model performance if its value exceeds 0.8 (Luo et al., 2017; Castellanos et al., 2019). TSS, calculated as sensitivity (correctly predicted presence points) plus specificity (correctly predicted absence or background points) minus one, indicates the prediction is reliable if its values exceeds 0.7. To reduce uncertainty associated with individual models, a weighted average ensemble model was created using only those models with an AUC value greater than 0.9, implying good predictive accuracy for final prediction (Marmion et al., 2009).

We first calculated the percent contribution of each predictor for the distribution of each species using the permutation method provided by the SSDM package. We then categorized the predictor variables into three groups as previously mentioned and assessed the contribution of each group to the distributions of pelagic and benthic species separately. Because the substrate variables were assumed to be temporally static, we did not analyze the contribution of these variables to the distribution of species. We employed Variation Partitioning Analysis (VPA) to partition the contributions from the physical oceanographic and physiology-based groups into individual and common contributions to the distribution patterns of species (Yang et al., 2023b). For the most important variable, water temperature, we fitted a response curve for each species using the GAM.

We predicted the present-day distribution of each species and its shift in the face of future climate change under three climate change scenarios in the 2040-2060 and 2080-2100. The prediction from the ensemble model is the occurrence probability of species with the value between 0 and 1, indicating the least and best habitat suitability. To better show the habitat distribution of each species, we converted continuous habitat suitability maps into binary format by maximizing sensitivity and specificity equality using the SSDM package (Liu et al., 2013). Species richness per grid cell was calculated by aggregating presence-absence maps of multiple species. Additionally, binary distribution maps for benthic and pelagic species were summed separately to assess species richness across different marine strata (Fitzpatrick et al., 2008). The standard deviation ellipse (SDE) is an effective spatial statistical method for accurately revealing the spatial distribution characteristics of geographic features. The center of the ellipse represents the relative position of distribution center, while the major axis direction indicates the expansion direction of distribution (Ren et al., 2022; Duman et al., 2023). We used SDE to characterize the distribution center and range of species under different climate change scenarios. To compare the distribution shifts among species with diverse ecological preferences, we calculated the proportional changes in their distribution areas across different climate scenarios. Additionally, we estimated the future spatial movement direction and distance of distribution center for each species using the SDMToolbox 2.0 package (Brown, 2014).

The model performances evaluated based on AUC and TSS for each individual model revealed variability among species ( Supplementary Tables 5 , 6 ). After filtering and weighting, the ensemble models for all the species exhibited improved performances compared with individual models, indicating by the overall increase of AUC and TSS values. The high AUC (above 0.9) and TSS (above 0.8) values implied the good predictive performance of the ensemble models for all the species ( Supplementary Table 4 ).

The percentage contributions of predictor variables provided by the SSDM model indicate that mean ocean temperature is the most significant factor in shaping the distribution patterns of benthic species, except for Pacific cod ( Figure 2 ). Depth makes the largest contribution to the distribution of pacific cod, followed by mean dissolved oxygen concentration and mean ocean temperature. Depth contributes the most to the potential distribution of pelagic species, followed by the mean ocean temperature. Altogether, mean ocean temperature plays a substantial role in the distributions of both benthic and pelagic species.

The results of the VPA indicate that 75% of the variation in the distribution of benthic species can be explained by the predictors from the physical oceanographic and physiology-based groups. The physiology-based group contributes the most, explaining 65% of the variation, while the physical oceanographic group explains 30%, and 20% of the contribution is shared by these two groups ( Figure 3A ). For pelagic species, 72% of the variation in their distribution can be explained by these two groups of predictors, with the physiology-based group again accounting for the majority, explaining 60% of the variation, while the physical oceanographic group accounts for 36%, and 24% of the contribution is shared by these two groups ( Figure 3B ).

The continuous and binary distribution of habitat suitability for each species under current climatic conditions is shown in Supplementary Figure 2 . The overlay of binary distributions over all species shows the distribution of species richness ( Figure 4 ). Areas with high species richness are concentrated along the coast of the Northwest Pacific Ocean, particularly in the shallow waters of the Japan Sea, the Bohai Sea, the Yellow Sea, and the East China Sea ( Figure 4 ).

Compared to the potential distribution areas at present day, areas of high numbers of species will continue to be located along the coast of the Northwest Pacific Ocean in the future. With the increase of greenhouse gas emission concentration, species richness per grid cell exhibits an expected trend of decrease at lower latitudes and increase at higher latitudes. The nearshore areas of the Yellow Sea, the East China Sea, and the Japan Sea will lose many species in the future. However, the coasts of the Okhotsk Sea and the Bohai Sea, as well as some offshore areas at higher latitudes, will attract more species to live in the future. Additionally, the spatial range of change in species richness will expand with the intensification of climate change ( Figure 5 ).

The results from the standard deviation ellipse analysis show that the main distribution areas of high species richness are located along the Northwestern Pacific coast, centered in the southern part of the Japan Sea, close to the Korea Strait. As climate change progresses, the distribution range gradually expands and the center moves northeastward ( Figure 6 ). The movement distance of the distribution center increases with time and intensity of climate change. It moves northeast by 104.06 km to the southeast area of Busan in the period 2040-2060, and then continues to move northeast by 168.48 km to the southern part of the Sea of Japan under the SSP5-8.5 scenario ( Supplementary Table 7 ).

From the species category perspective, our projections indicate a general decline in suitable habitat areas for crustaceans, including Acetes japonicus and Portunus pelagicus. This decline is more pronounced under higher emission scenarios (SSP5-8.5), where both species exhibit substantial contractions in their habitats by 2080-2100. In contrast, changes in fish habitat areas under future climate change depend on their ecological preferences. All the tropical and temperate fishes in our study, including Harpadon nehereus, Trachurus japonicus, and Engraulis japonicus will experience habitat expansion compared to their current habitat. Two subtropical species show opposite trends in habitat change, with an expansion of Scomber japonicus and a contraction of Larimichthys polyactis in the face of future climate change. Two boreal species also show opposite change trends of habitat, with a great expansion of suitable habitats for Clupea pallasii and a slight contraction of suitable habitats for Gadus macrocephalus. Considering the living water strata of species, all four pelagic species are projected to expand their habitats under three climate change scenarios for both periods in the future. In contrast, all five benthic species, except for Harpadon nehereus, will contract their habitats under three climate change scenarios for period 2080-2100 ( Table 1 ; Supplementary Table 8 ).

The distribution centers of species are linked to their respective climatic zones and typically shift towards higher latitudes as climate change progresses, with migration distance increasing with the intensity of climate change. Under present climatic conditions, Gadus macrocephalus, adapted to boreal climates, exhibits the northernmost distribution, followed by Clupea pallasii. In contrast, the tropical species Harpadon nehereus shows the southernmost distribution. The calculated results of movement distances of the distribution center for each species indicate that crustaceans typically migrate shorter distances than most fish species in response to climate change. Among fish, pelagic species typically exhibit greater migration distances than benthic species. Clupea pallasii shows the longest average migration distance, approximately 431.52 km, followed by Scomber japonicus with an average of 326.37 km. In contrast, the benthic species Harpadon nehereus and Gadus macrocephalus demonstrate shorter average movement distances, at 82.84 km and 287.47 km, respectively. Considering thermal preferences, fish inhabiting boreal and temperate climate zones with higher latitudes demonstrate greater migration distances than those inhabiting at tropical and subtropical zones with lower latitudes, implying that their habitats are more inclined to be impacted by climate change. Notably, the subtropical species Scomber japonicus exhibits a migration distance of 1115.95 km in the periods 2080-2100 under the SSP 5-8.5 scenario that far exceeds other tropical and subtropical fishes ( Figure 7 ; Table 3 ).

Pelagic species seem to be more resilient to changing climate than benthic species, as evidenced by their better ability to colonize new areas in the future. The stronger swimming or dispersal ability of pelagic species facilitates these species to occupy more regions, resulting in a broader range of spatial distribution and environmental niche. Studies have found that climate change facilitates habitat expansions of environmentally tolerant species (Thomas et al., 2016). For instance, our study showed that Scomber japonicus can live within a temperature range of 5-25°C ( Supplementary Figure 4 ), and this broader adaptive range facilitates this species to expand its habitat in the future (Han et al., 2023).

Benthic crustacean species commonly have limited mobility and do not move or move only a short distance in a lifetime. The narrow ranges of spatial distribution and environmental niche may result in vulnerability to future climate change (Green et al., 2014; Cloyed et al., 2021; Dubos et al., 2022). Our study indicates a decrease in suitable habitat for Acetes japonicus and Portunus pelagicus and the reduced extent increases with environmental change intensity. The largest habitat reduction was found in the period 2080-2100 under the SSP5-8.5 scenario. Compared to benthic crustaceans, benthic fish usually possess superior swimming abilities, facilitating them to adapt to a wider range of environmental conditions. However, we found a habitat contraction for Larimichthys polyactis and Gadus macrocephalus, consistent with published studies (Li et al., 2022; Yang et al., 2023a), and a habitat expansion for Harpadon nehereus. This contrasting result implies that climate-induced changes in the benthic environment might mostly compress the distribution ranges of species, especially for cold water species, but can also generate more suitable habitats, especially for warm water species. Wang et al. (2021) predicted a habitat reduction of Harpadon nehereus in the coast of China due to climate change. We used a greater number of predictor factors with updated environmental data to predict a habitat expansion at a larger spatial scale. The differences in data and spatial scale may be the reasons for generating this inconsistent result. Our results indicate that for the cold-water species, Pacific cod (Gadus macrocephalus), habitat area may experience a slight increase under the SSP2-4.5 scenario during the 2040-2060 period compared to other climate scenarios. This increase is likely driven by moderate warming in high-latitude regions, where previously colder areas become more conducive to species distribution ( Supplementary Figure 3 ). Moderate warming extends the range of favorable environmental conditions in these regions, thereby enhancing habitat suitability. Conversely, under SSP1-2.6, which represents lower warming, the temperature increase may be insufficient to significantly expand suitable habitats in colder regions, resulting in a more limited range expansion. For Portunus pelagicus, our results indicate that the SSP2-4.5 and SSP5-8.5 scenarios, compared to SSP1-2.6, may result in the acquisition of more suitable habitats at higher latitudes ( Supplementary Figure 3 ). This is because the Portunus pelagicus, a demersal species with an optimal temperature of 26°C ( Supplementary Figure 4 ; Liao and Li, 2002), benefits from slower rates of environmental change in the demersal layer, which is less immediately affected by climate change (Heuzé et al., 2015). Climate change could gradually shift the species’ suitable habitats to higher latitudes as it seeks optimal environmental conditions.

In terms of climatic zone, we expect that tropical fishes may be more resilient to climate change than fishes in other climate zones because large amounts of suitable habitats may arise in temperate or boreal zones (Nakamura et al., 2013; Osland et al., 2021; Haslam et al., 2023). We predicted habitat expansions and northward movements of the distribution center of tropical species, Trachurus japonicus and Scomber japonicus, in the face of climate change. However, this result may change when the spatial extent of the study area changes. For instance, Xiong et al. (2024) estimated the impact of climate change on the distribution of Trachurus japonicus in the Northern South China Sea and found the distribution of this species expanded northward and contracted southward. The rate of habitat contraction surpassed that of expansion, with an estimated decrease in suitable habitats by the end of the 21st century. When we extended the study area, especially to the north, it is possible to find that the rate of habitat expansion surpasses that of habitat contraction.

In terms of fish size, studies indicated that small pelagic fishes exhibit stronger adaptability to climate change than big fishes (Lefort et al., 2015). Similar to previous findings (Montero-Serra et al., 2015), our results indicate that climate change will drive the small pelagic fish Engraulis japonicus to expand its distribution in the period 2080-2100 under the SSP2-4.5 and SSP5-8.5 scenarios. A long-term change of environmental conditions can gradually motivate genetic variation to adapt to environmental change stresses (Geerts et al., 2015; Scheffers et al., 2016; Krajmerová et al., 2017). Small fishes usually have a short lifespan and reproductive cycle and thus a short time interval in the alternation of generations, which may speed up the genetic variation to adapt to climate change (Munday et al., 2013). Furthermore, Benedetti et al. (2021) predicted that the abundance of diatoms, dinoflagellates, haptophytes, and other plankton will increase in temperate zones, leading to a rise in the biomass of small fish species feeding on them. Climate change is also expected to drive changes in plankton community composition, with a shift towards smaller species at higher latitudes (Henson et al., 2021). This shift could significantly impact the nutritional value and availability of specific prey species, potentially affecting the long-term sustainability of species reliant on them. The strong environmental tolerance and richer food sources induced by climate change are among the reasons for the habitat expansion of small fishes, such as Engraulis japonicas.

We separated predictor variables into three groups, and physiology-based group showed dominant importance in affecting the distribution patterns of species in the Northwest Pacific Ocean. This implies that climate-induced changes in the physiology-based factors may result in physiological maladjustment and then prompt adaptive evolution of species, such as alterations in metabolism, growth patterns, and reproductive strategies (Bates et al., 2014). Though this adaptive evolution differs among species, it may alleviate climate change impact on species.

Our study demonstrated that the distribution centers of all nine species, whether expanding or contracting, will shift northward as a result of climate change. This redistribution may diminish the effectiveness of existing spatial fisheries management strategies. For example, our predictions suggested that Larimichthys polyactis, an economically important species in China, will lose substantial suitable habitats in the East China Sea by 2080-2100 under the SSP5-8.5 scenario. This suggests that current management practices, such as the establishment of national aquatic germplasm resource conservation areas in 2008 to protect spawning grounds for hairtail (Trichiurus haumela) and migratory routes of small yellow croaker (Larimichthys polyactis), will become less effective by the end of the 21st century.

Furthermore, these distribution shifts of economically significant species may make transboundary fisheries face more serious challenges because these fisheries are prone to over exploitation even if climate change effects on stocks were taken into account (Liu and Molina, 2021). The inconsistent management policies between nations may exacerbate these challenges. Our predictions suggested that Pacific cod (Gadus macrocephalus) habitats will shrink in the China Seas and Sea of Japan while increasing in the Sea of Okhotsk. This distribution shift is unfavorable to China, Japan, and Korea, while benefiting Russia. Similar scenarios will occur for other species in various regions (Astthorsson et al., 2012; Ojea et al., 2020). These new distribution patterns induced by climate change underscore the necessity for strengthening international cooperation to manage transboundary stocks. Regional Fisheries Management Organizations (RFMOs), such as the South Pacific Regional Fisheries Management Organization (for jack mackerel), will face growing pressure to adapt their frameworks to prevent overexploitation and balance stakeholder interests on shared stocks.

Spatial fisheries management tools are critical for the sustainability of stocks. For instance, Liu et al. (2018) applied these tools to address stock migration on groundfish fisheries, offering a potential model for adapting to climate-induced distribution shifts. While long-term fishery resource surveys are critical for designing fishery management strategies, predictive models can predict future distribution of species under different climate scenarios, allowing managers to explore and plan for possible changes in biomass or distribution, rather than being caught on the back foot once it has already happened. Therefore, our models can be as part of a toolkit that fisheries managers can use to ensure sustainability of fish stocks into the future.

Though our study is significant in informing climate-smart fishery management strategies in the Northwest Pacific, there are several limitations that need to be addressed in the future work. First, although our modeling results indicate that species with different ecological preferences exhibit varying distributional responses to climate change, the intrinsic assumption of predictive models is that the niche space of each species does not change over time (Guisan et al., 2017; Smith et al., 2019). Ignoring intraspecific climate adaptation and population response variability to climate change may diminish prediction accuracy (Hällfors et al., 2016). Mechanistic models that incorporate population genomics are suggested to offer more reliable estimates than conventional statistical models. Overlooking adaptive evolution may lead to an overestimation of climate impacts on species’ habitats (Robinson et al., 2011; Razgour et al., 2019). Second, our predictions assumed that all species have unrestricted dispersal capabilities. This over-optimistic assumption implied that each species can migrate to all suitable areas under climate change. However, ecological and man-made barriers and species mobility may limit some species from reaching suitable habitats (Tamario et al., 2019; Liang et al., 2021). Therefore, this assumption undoubtedly overestimated species’ suitability to climate change. Third, due to data limitations, our study considered only environmental influences, excluding biological factors and interactions among them such as predation, competition, parasitism, mutualism, and facilitation. These relationships can interact in complex ways with climate change, influencing species’ spatial distribution patterns (Tylianakis et al., 2008; Braschler and Hill, 2007). For example, Estes et al. (2011) suggest that large apex predators can significantly limit the distribution and range of prey species in marine ecosystems. The exclusion of such interactions from our models may lead to an overestimation of available habitats for prey species, as predatory pressures are not considered. Fourth, most economically important species inhabit coastal continental shelf waters, where commonly face heavy fishing pressure and pollutions from human activities (Eigaard et al., 2017; Palummo et al., 2023). Ignoring fishing pressure and other disturbances by human activities will undoubtly overestimate the distribution range of species.