The research area of this study was set along the coast of subtropical to temperate regions of Japan between 24.0°N and 41.6°N and 122.9°E and 142.2°E (Figure 1). From the 4,500 fish species that inhabit in Japanese waters (Motomura, 2020), we selected six fish species that represent different ecosystem services to humans and that have sufficient occurrence data (n > 10) for the estimation of recent and future distribution by the maximum entropy model (MaxEnt) (Wisz et al., 2008; Ishihama, 2017). We selected Aluterus scriptus and Scarus ovifrons as toxic harmful fish, Kyphosus bigibbus and Siganus fuscescens as herbivorous fish, and Amphiprion frenatus and Chaetodon auriga as tropical reef fish (Table 1). These toxic harmful fish contain a kind of palytoxin in their internal organs. These species could be mistaken for Stephanolepis cirrhifer and Chlorurus microrhinos because of their similar shapes. Food poisoning can cause death by accidental ingestion, and poisoning has been reported in more than 100 people between 1953 and 2020 (Ministry of Health, Labor and Welfare Japan, 2021). Distributional changes in these species increase the risk in unfamiliar regions (Ministry of Health, Labour and Welfare Japan, 2021; Taniyama et al., 2003). In addition, expanding of herbivorous fish causes barren grounds of seagrass/seaweed beds due to their increased grazing pressure under warm water temperatures (Vergés et al., 2016; Kumagai et al., 2018; Fisheries Agency Japan, 2021). On the other hand, these tropical reef fish are popular for marine leisure and aquarium exhibition purposes.

The distribution data of six targeted fish species consisted of a total of 871 records collected and stored in open source databases using the GBIF (2021) (https://www.gbif.org), OBIS (2021) (https://obis.org/), BISMaL (https://www.godac.jamstec.go.jp/bismal/j/), Science Museum Net (S-Net; http://science-net.kahaku.go.jp/), and the Fish Image Database (2021) (FishPix; http://fishpix.kahaku.go.jp/fishimage/). For the FishPix data, we obtained the coordinate information with third decimal coordinate accuracy using the site name from the list of diving spots for the sighting information data only with the site name (Marine Planning, 2001; Professional Association of Diving Instructors Diving Spot Information, 2021). We corrected spatial biases to one or two occurrence data per 0.5° grid (Boria et al., 2014; Fourcade et al., 2014). A total of 156 of these data for the six indicator fish species were recorded between 1962 and 2021 (Table 1). We excluded young or juvenile fish occurrence data because of the possibility of abortive migration that appears from summer to autumn and disappears during winter in warm current areas (Shimizu et al., 2009; Booth et al., 2011), although some databases do not provide life stage status. Thus, this dataset includes abortive migration fish due to insufficient life stage data. For quality control of the dataset, we also excluded low accuracy occurrence data because it was lower than the third decimal coordinate accuracy.

Species distribution models (SDMs) identify the relationships between the presence of a species and environmental variables. Many algorithms exist for SDMs, among which MaxEnt is the most commonly used because it is robust against georeferencing errors using presence-only records (Phillips et al., 2006; Graham et al., 2008; Elith et al., 2011). It outperforms most other algorithms, such as generalized linear models, generalized additive models, and random forests, and it is particularly appropriate for marine species (Elith et al., 2006; Ready et al., 2010; Aguirre-Gutiérrez et al., 2013). MaxEnt also performs well in estimating potential range shifts for a species due to climate change (Hijmans and Graham, 2006). Therefore, MaxEnt was selected for this study. All MaxEnt models were run using the default settings (version 3.3.3 k; Phillips et al., 2006) and replicated ten times using cross-validation (Marsland, 2009). The entire data were randomly divided into ten parts: 10% of the data were used for verification, and the remaining 90% were used for learning. All the divided parts were used for verification during the ten runs of the model to evaluate the model results. The model accuracy was examined using the area under the curve (AUC) (Fielding and Bell, 1997), for which values >0.7 are commonly accepted (Swets, 1988; Raes and ter Steege, 2007). Predicted logistic values of each grid were converted to presence/absence values using a ten percentile training presence logistic threshold. This threshold was calculated based on the omission errors, which allow a 10% false presence of the training data. Logistic values lower than this threshold were classified as absence, and higher probabilities were classified as presence (Escalante et al., 2013).

We created a database of five environmental variables at 1 km resolution. We chose the following variables: 1) coldest month sea surface temperature (SST), 2) depth, 3) slope, 4) seagrass/seaweed bed area, and 5) coral reef area (Table 2). Cold month sea surface temperature was chosen because cold temperature restricts the distribution range of tropical fish. Depth was used because the variable well explained coral reef fish and temperate reef fish habitat. Slope represents topographic complexity that affects fish habitat preference, and seagrass/seaweed bed area well-explained the occurrence of nearshore temperate reef fishes (Pittman et al., 2007; Young and Carr, 2015; Leriorato et al., 2021). The coral reef area was used as a proxy for tropical reef fish habitat. We set the resolution of all variables to a 1 km grid by interpolation with the inverse distance weighting method (Bartier and Keller, 1996). Among the various seasonal SST data, we used the coldest month water temperature because it has been discussed that the distributional changes of tropical fish are most sensitive to the lowest temperature (Leriorato et al., 2021). This is because an abortive migration of larvae occurs from summer to autumn and disappears in winter due to low temperatures in the warm current area (Nakazono, 2002; Shimizu et al., 2009). The coldest month SST data were obtained from Modis-Terra with an approximately 4-km horizontal resolution and averaged for the period of 2000-2018 to avoid the effect of the Pacific decadal oscillation (Mantua and Hare, 2002). The seafloor structure, such as depth and topography, is one of the important factors to explain a suitable habitat for fish (Pittman et al., 2007; Costa et al., 2014; Young and Carr, 2015). For the depth data, we used Digital Bathymetric M7000 which is the most accurate data around Japan, and averaged depth data in each 1-km grid. We then calculated the slope from the difference in depth and distance between the grid and eight neighbors (Tan θ = difference of depth/distance) using the ESRI ArcGIS Spatial Analyst slope function. The complex structure of coral reefs provides a habitat for many tropical reef fish species. On the other hand, the seaweed and seagrass provide food resources for many fish species. Therefore, occurrence of two target tropical reef fish species and herbivorous fish is likely related to the presence and absence coral reefs and seaweed/seagrass beds. For these reasons, we included the coral reef area and seaweed/seagrass bed area in each 1-km grid as the environmental variables. The coral reef dataset included those in the temperate region along the coast of the Pacific Ocean affected by the warm Kuroshio Current. The seaweed/seagrass dataset was based on satellite image analyses, which evaluates the area of large seaweeds (such as sargassum and kelp) and seagrass.

The future ocean regional projection (FORP) dataset was produced by high-resolution regional ocean model simulations with ensemble atmospheric forcing from Coupled Model Intercomparison Project phase 5 (CMIP5; Taylor et al., 2012) models and with a future emission scenario based on the representative concentration pathways (RCPs; Moss et al., 2010). These were developed by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) under the Social Implementation Program on Climate Change Adaptation Technology (SI-CAT) (Nishikawa et al., 2021). The Meteorological Research Institute Community Ocean Model version 4 (MRI.COMv4; Tsujino et al., 2017) was used for the regional ocean models. FORP-JPN02 version2 is a historical and future ocean projection dataset surrounding Japan with a horizontal resolution of approximately 2 km, produced by the dynamical downscaling simulation. This dataset precisely expresses coastal topography and ocean currents to apply to coastal areas. We used the averaged results from four climate models from FORP-JPN02 version2 (GFDL-ESM2M, MRI-CGCM3, MIROC5, and IPSL-CM5A-MR) to improve the performance and decrease the uncertainties of a forecast based on the fundamental assumption that climate model results generally have biases and errors (Tebaldi and Knutti, 2007). We used climate model results for 2046–2055 under RCP8.5 and during 2091-2100 under RCP2.6 and 8.5 simulations to examine the effects of ocean warming on potential habitats. RCP2.6 is the scenario that aims to keep global warming likely below 2°C above pre-industrial temperatures, and RCP8.5 scenario is the worst-case scenario assuming greenhouse gas emissions continue to rise throughout the 21st century (IPCC, 2014). Other RCP simulations have not yet been provided for this dataset. The averaged results from four climate models (for the years 1996-2005) in FORP-JPN02 version2 were 0.91 ± 1.06°C lower than the satellite observation data for 2000-2018. The biases between the observed and the model results were corrected by adding the difference of the model results to the observed climatology using the method of Yara et al. (2011). We added the difference values between the model results from 1996–2005 and 2046–2055 under the RCP8.5 simulations to the satellite observation data of SST and between 1996–2005 and 2091–2100 to the observed SST under the RCP2.6 and RCP8.5 simulations. SST data from 2000–2018 were substituted with RCP climate scenario data for the projection models, and suitability was recalculated using MaxEnt for each species. The shift in the latitudinal distribution of each target fish species was compared between the Pacific coast and the Sea of Japan. The total grid number within each degree of latitude (e.g., 29.00°N to 29.99°N) was calculated by summing the number of grids with a species’ presence. Even though it is not directly influenced by the Kuroshio Current, the coastal area of the East China Sea was classified as the Pacific coast in this study for graphical representation. Water depth, seagrass/seaweed areas, and coral reef areas will change in response to increase in water temperature. Although the sea level is anticipated to rise in accordance with global warming (i.e., Oppenheimer et al., 2019), a precise projection of sea level rise is difficult by using models with a spatial resolution of ∼1 km, as used in this study. Furthermore, to predict precisely the future distribution of coral reefs and seaweed beds, various factors such as water temperature, ocean acidification, the rise in sea level, and wave height, are needed to be considered (Yara et al., 2012; Takao et al., 2015; Morim et al., 2019). Therefore, it is difficult to predict their changes in this study. We, therefore, used the same water depth, seagrass/seaweed, and coral reef area data of the 2000–2018 estimation for future predictions.

The predicted distribution of six fish species in 2000–2018 varied from 4,006 to 16,430 grids in the models (Table 3). Among them, S. ovifrons, C. auriga, and A. frenatus had the same northern limit, ranging from 35°N along the Pacific Ocean coast. In addition, semi-closed seas such as Seto Inland Sea, Ise-Mikawa Bay, and Tokyo Bay were not suitable for these species on the Pacific Ocean side (Figures 2–7, Supplementary Figure S2). A. scriptus had the widest distribution ranging from 38°N, along the Sea of Japan. The second most abundant species K. bigibbus and S. fuscescens were estimated to occur to 36–37°N along the Sea of Japan. On the other hand, A. frenatus, S. ovifrons, and C. auriga were limited to 34°N along the Sea of Japan (Figures 2–7, Supplementary Figure S4).

The distribution of the six fish species was well predicted by the combination of five environmental variables with a mean AUC of 0.82 to 0.93 in the models (Table 4). The SSTs accounted for nearly 50% to more than 80% of the SDMs. In contrast, the slope was the second most influential factor for S. ovifrons and K. bigibbus. Depth was the second most influential factor for A. scriptus and S. fuscescens that inhabit shallower depths. Seagrass/seaweed beds contributed 7–8% of the estimated distribution for S. ovifrons and S. fuscescens.

All species were predicted to expand in distribution under the two different future carbon emission scenarios (Table 3). The model with RCP8.5 by the 2091–2100 forecasted wide expansion to northeastwards (Figures 2–7, Supplementary Figure S4). The predicted expansion differed between the Pacific coast and the Sea of Japan.

All species were predicted to expand in distribution under the two different future carbon emission scenarios (Table 3). The model with RCP8.5 2091–2100 forecasted wide expansion to northeastwards (Figures 2–7, Supplementary Figure S4). The predicted expansion was different between the Pacific coast and the Sea of Japan.

The predictions showed significant differences in the northern shift among species on the Pacific Ocean side. For C. auriga, A. frenatus, and S. ovifrons, the model RCP2.6 and 8.5 predicted a slight northern shift along the coast of the Pacific Ocean by 2091–2100. Similarly, for A. scriptus, K. bigibbus, and S. fuscescens, the model RCP2.6 predicted a slight northern shift along the coast of the Pacific Ocean by 2091-2100, and the model RCP8.5, by 2046–2055. In contrast, the model RCP8.5 predicted a 500–600 km northern shift along the coast by 2091–2100 (Figures 2–7, Supplementary Figure S4).

Some species were predicted to expand in semi-closed seas located on the Pacific Ocean side (Supplementary Figure S2). For S. fuscescens, the model RCP8.5 was predicted to expand to Seto Inland Sea, Ise–Mikawa Bay, and Tokyo Bay by 2091–2100. In contrast, for A. scriptus, K. bigibbus, A. frenatus, and C. auriga the model RCP8.5 predicted to expand to only Seto-Inland Sea by 2091–2100. In contrast, for S. ovifrons, the model predicted expansion only along the coast of the Pacific Ocean.

The predictions showed significant differences in the northern shift in the Sea of Japan compared to the Pacific Ocean. For all the target species, the model RCP2.6 predicted a 100–300 km shift north by 2091–2100 along the coast of the Sea of Japan, and also the model RCP8.5, by 2046–2055. However, the model RCP8.5 predicted 200–400 km by 2091–2100. Among the six species, only a 100–200 km shift north was predicted for A. frenatus for both RCP2.6 and 8.5 scenarios by 2091–2100.

The coldest month SSTs were selected as key environmental factors for the models in our study, which agrees with previous studies on S. fuscescens distributional change in Japan (Takao et al., 2015). Water temperature is an important environmental factor for ectotherms. Body temperature is affected by environmental water temperature in most fish species and controls the survival, feeding, growth, reproduction, and geographical distribution of organisms (Angilletta et al., 2006). Physiologically, fish survive in a particular temperature range, and natural geographical distribution is related to minimum water temperature (O'Connor et al., 2007; Payne et al., 2016). An increase in the water temperature affects the starvation risk associated with a high metabolic rate (Houde, 1989; Munday et al., 2012). Food consumption of K. bigibbus was the highest at 30°C, declined to less than 16°C, and stopped at less than 12°C in the aquarium experiment (Noda et al., 2016). The feeding activity of S. fuscescens becomes most active at 26–29°C, decreases at temperatures lower than 20°C, and stops at < 15°C (Fisheries Agency Japan, 2021). In addition, A. scriptus survives at temperatures higher than 14°C (Tanase et al., 1992), and C. auriga survives at 15–33°C (Tsuchida et al., 2018). It is likely that the distribution of fish studied here is mostly sensitive to changes in temperature with the ongoing climate change.

Seafloor structures, such as depth and topography, are important factors in explaining suitable habitats for many fish species (Pittman et al., 2007; Costa et al., 2014; Young and Carr, 2015). In this study, depth and slope contributed to the model results and correlated with the lower depth and steeper slope for the target species (Supplementary Figure S3). Adults of A. scriptus are seen along deep coastal slopes or outer reef drop-offs at a depth of around 20 m, and they are usually captured between 3 and 20 m (Kuiter and Tonozuka, 2001; Luiz Jr et al., 2008; Izzo et al., 2010). Their habitat preference was consistent with depth and slope contribution of the model. Some fish change their habitat in a day because of feeding and hiding. For example, S. ovifrons spends most of their time in shallow water lower than 10 m in the daytime. However, S. ovifrons sleeps at 10–15 m during nighttime (Gomi et al., 2021). Biotelemetry attached to K. bigibbus showed that the most time was spent at a depth of 10–15 m near the seafloor with migration to shallow water to search for seaweeds in western Japan (Yamaguchi et al., 2006). These behaviors show a habitat preference for steep-slope along rocky shores.

Coral reefs provide habitats for tropical reef fish. Thus, we hypothesized that coral reefs contribute to these models especially for A. frenatus and C. auriga. However, coral reef contributions were regarded to be low for all species by the models. Coral reef fish can survive not only in coral reefs but also in rocky reefs, given the availability of food resources. C. auriga feeds on benthic invertebrates, and A. frenatus inhabits specific anemones and does not inhabit corals (Hattori, 2011; Nakamura et al., 2013).

The presence of seagrass/seaweed was selected as an important factor for S. ovifrons and S. fuscescens, which have different food preferences. Scarus ovifrons is a scraper that feeds on small epiphytic algae on rocky substrata and dead corals (Steneck et al., 2017; Gomi et al., 2021). The contribution of seagrass/seaweed beds to the S. ovifrons model can be attributed to the fact that the macroalgae vegetate on rocky reefs and the rocky reefs form the habitats of S. ovifrons. In contrast, S. fuscescens is a browser that prefers to graze seaweeds and seagrasses (Noda et al., 2011). For these species, the seagrass/seaweed was positively correlated with the probability of occurrence. In contrast, the seagrass/seaweed did not contribute to K. bigibbus, the main diet of which is seaweed. This species was not found in seaweed beds during the daytime. However, their groups were found around structures such as breakwaters and rocky areas (Kadota et al., 2017). Their specific behavior can cause low contribution of seaweed beds for this species.

The high contribution of minimum water temperature for all the target species was caused by the wide range of evaluations stretching from the southern to northern part of Japan. Additionally, the depth and slope were contributing factors for all the species because none of them were found in tidal flats, gentle slopes, and deep waters. However, these results may have sampling bias because some areas, such as steep slope or deep water, are difficult to survey. For a more strict evaluation, data on true-absence of the species in those areas need to be collected.

The AUC is the probability that a randomly chosen datum of presence is ranked above a datum of pseudo-absence chosen randomly from the study area (Phillips et al., 2006). Basically, modeling on the basis of presence only renders evaluation of the model accuracy difficulty, owing to the lack of true-absence data. Machine learning models based on presence only tend to overfit the model accuracy (Warren and Seifert, 2011; Merow et al., 2013). The model results include accidental noise in the occurrence data, which decreases the generality of the model results and its predictive ability (Wenger and Olden, 2012). If there is a sufficient amount of less-biased occurrence data, it would be possible to reduce the overfit by using the different dataset for the evaluation. Moreover, low-resolution environmental variables also tend to show high accuracy because the relationship between occurrence data and environmental variables becomes unclear owing to the averaged data in a large grid (Randin et al., 2009). Although a 1-km resolution is high for large-scale modeling, this study has a possibility of overestimating the AUC.

Ongoing water temperature rise with global climate change is predicted to shift fish distributions toward higher latitudes, as revealed in this study and others (Takao et al., 2015; García Molinos et al., 2016). In fact, shifts in the global distribution of fish species have already been observed in the past (Nakamura et al., 2013; Vergés et al., 2014; Wernberg et al., 2016; Kumagai et al., 2018). These fish have overcome abortive migration and will expand their distributional range according to future increases in water temperature.

The expansion of the distributional range was different between the Pacific Ocean and the Sea of Japan side. This difference can explain in relation to ocean current systems around Japan. The indicator fish species can expand the distribution range on the Pacific Ocean side, which is influenced by a warm Kuroshio Current (Figure 1). However, the distribution shift further north was limited due to the Oyashio Current that flows from the northeast to the southwest, meeting Kuroshio at around 35–36°N. In contrast, tropicalization expands to higher latitudes in the Sea of Japan side influenced by the Tsushima Current, which diverges from the Kuroshio Current around 30°N and flows along the coast up to 40°N or a higher latitude. Kumagai et al. (2018) also provided evidence of the interaction between water temperature increase and ocean currents for the dispersal of marine organisms using the 60-year records in this region. Although we attempted to minimize the biases and errors by using the ensemble model, the results of future projection models still include uncertainty. The relative increase in water temperature of the ensemble model for the years 2091-2100 was 1.2 ± 0.4°C for RCP2.6 and 3.6 ± 0.8°C for RCP8.5 (Supplementary Table S1).

Future projections of S. fuscescens were different from a previous study by Takao et al. (2015) who predicted a narrower expansion along the coast of the Sea of Japan by 2100 using similar scenarios. This discrepancy is ascribed to the use of different ocean models for forecasting future ocean conditions and the use of SST in different seasons. Because our study used higher resolution ocean model outputs (FORP-JPN02 version2) which also includes information on ocean current and topography and because we used the coldest water temperature which is more critical for determining the northern limit of tropical fish distribution, the results of our study provide a more realistic prediction of future shifts in the indicator fish species with multiple climate scenarios.

The models with RCP8.5 by the years 2091-2100 predict that there will be an expansion to semi-closed seas such as Seto Inland Sea, Ise-Mikawa Bay, and Tokyo Bay that are not directly affected by warm currents. Water temperature of these seas will increase at a slower rate than that in the coastal areas. Water temperature greatly contributed to the model for S. fuscescens which is widely distributed because of its tolerance to low temperatures (Figure 5, Supplementary Figure S1–2, Supplementary Figure S2, Supplementary Figure S3–2). Thus, tolerance to low temperatures is one of the factors for the expansion of S. fuscescens to semi-closed seas. Architects, Regional Planners and Associates et al. (2019) reported future predictions for S. fuscescens in the Seto Inland Sea. This report expanded the area of species range shifts by 90–100% by 2100 using the MRI-CGCM climate model (RCP8.5). The report applied 10.4–11.1°C as a lower limiting water temperature threshold which estimated survival of S. fuscescens in aquarium conditions (Ueta and Tanada, 2018, Supplementary Figure S1). Their study evaluated only the use of SST and applied the harshest threshold. There was natural population dispersal northward to cooler waters from breeding population in the south. These dispersal individuals have difficulty surviving in the winter due to low temperature (Lenanton et al., 2017). But in the summer, they overgraze the seagrass/seaweed. Thus, physiological performance in the aquarium often did not explain well the distribution of species in nature (Payne et al., 2016). Because of these reasons, the range shift was largely different from our model results.

The barren ground sites that were caused by herbivorous fish shifted northward from the 1970s to the 2000s and changed coastal ecosystem services in the western part of Japan (Kumagai et al., 2018). In addition, increases in water temperature enable these fish activities throughout the year and affect coastal ecosystem services. Suitable habitats for coral reefs and tropical seaweed/seagrass beds may also shift northward in the future. However, temperate seaweed species disappear as the water temperature rises (Yara et al., 2011; Tanaka et al., 2012; Takao et al., 2015). However, the dispersal velocity of algae is slow compared to that of other taxa (Poloczanska et al., 2013). If the northern shift velocity is faster than the dispersal velocity, it may cause the disappearance of suitable habitats.

We only evaluated six fish species (two for each of the three functional groups) as indicators due to the limitations of collecting adequate occurrence data. More fish in the same functional groups occur in the same distributional range, for example, Pterois lunulata and Hapalochlaena fasciata as toxic harmful species, Calotomus japonicus and Prionurus scalprum as herbivorous fish, and Chrysiptera cyanea and Nemateleotris magnifica as coral reef fish that are attractive to the tourists. These species are also expected to expand their distributional range northward, similar to our indicator species, affecting various ecosystem services for the human society. A precise prediction of future forecasts of these indicator species, derived from this study, is thus useful for preparing adaptation measures against changes in ecosystem services due to distributional changes of fish.

The region predicted as the suitable habitat for toxic harmful fish in the future would experience more cases of human food poisoning caused by consuming these fish and more injuries by accidental contact during marine leisure. To decrease such incidents, local governments may need to raise public awareness with fishers to prevent these species from being sold in the marketplace and with managers of beach and sea bathers to avoid these toxic fish that were previously unknown on their beaches. To mitigate the expansion of herbivorous fish, local governments and fishery communities in the western regions of Japan need to prepare for the effective protection of seagrass/seaweed beds. Currently, effective methods to remove these detrimental herbivorous fish species are not well developed. However, some efforts such as the exclusion of fish by setting bait cages in seagrass/seaweed beds and selective harvesting of these fish may be considered as management measures. Extraction from fisheries catch and promotion of the consumption of these detrimental herbivorous fish in the marketplace can reduce their population and negative effects (Fisheries Agency Japan, 2021). Suitable site selection for seaweed aquaculture can assist local fisheries (Liu et al., 2013). On the other hand, an increase in the abundance and species richness of tropical reef fish may attract tourism, such as snorkeling and scuba diving, which can benefit local economies. In addition, the monetary value of fisheries production and the species richness of aquarium fish may increase in the temperate regions in the future, and ecosystem services from the reef fish may change (Biodiversity Center of Japan, 2015; Sato et al., 2021). Changes in the distribution of coral-dependent tropical fishes, such as obligate coral dwellers and corallivores, cannot be estimated using this model. This is because changes in the coral reef distribution were not considered in this model. Therefore, the estimated changes in the distribution of the two tropical fish species do not indicate changes in the future distribution of all coral reef fishes.