Influences of Global Warming on the Larval Survival and Transport of Snow Crab ( Chionoecetes opilio ) in the Sea of Japan

The snow crab (Chionoecetes opilio) sustains an important bottom trawling fishery in the Sea of Japan. Its response to global warming is attracting the attention of the public. Using a transport and survival model for crab larvae in the Sea of Japan, we examined the spatial-temporal variations of crab spawning and larval settlement in the past (mid-20th century), present (early 21st century), and future (midand late 21st century) under the low and high radiative forcing scenarios. It was found that the variations in spawning differed between the regions south of and north of 41.5◦ N, on both seasonal and long-term scales. Larval settlement in the Sea of Japan was projected to increase in the future, which is mainly attributed to a reduction in mortality due to the low water temperature. Moreover, the aggregating location of the settled megalopae will likely shift northward, with increasing settlement off Hokkaido Island. With additional sensitivity experiments, we confirmed that the change in water temperature has a stronger impact on larval settlement than that in the current field. The change in water temperature controlled both the amount and distribution of crab larval settlement, while a change in current field only affected the distribution to some extent.


Introduction
The snow crab Chionoecetes opilio (C.opilio) is distributed in deep cold waters in the northwest Atlantic Ocean and north Pacific Ocean, including the Sea of Japan.It is an important commercial fishery species, particularly in Japan, Alaska, and Atlantic Canada [1][2][3][4][5].Snow crabs hatch between early February and late April, with a peak in March.After three larval stages (first, second zoea, and megalopa), larvae settle on the seabed and metamorphose into the first stage benthic crab in June [6].Female crabs molt eleven times in 7-8 years before they are mature and reproductive.Male crabs complete nine molts and mature one year earlier than females, but they do not mate until they undergo the terminal molt between the eleventh and sixteenth molt [7,8].As the most important species in the offshore bottom trawling fishery in Japan, the catches of snow crab declined dramatically during the mid-to late 20th century, mainly due to overfishing [9].However, changes in marine environmental conditions, due to climate change, have also contributed to the decline [10,11].
Historical observational data have shown a trend in global warming since the second half of the 19th century.Global annual mean surface air temperature has increased at a rate of 0.72 • C/century and sea surface temperature (SST) increased at a rate of 0.53 • C/century [12].In the projection of the future climate, a set of four Representative Concentration Pathways (RCPs) was developed, representing different socio-economic and emission scenarios [13], including one mitigation scenario leading to a very low forcing level (RCP2.6),two medium stabilization scenarios (RCP4.5/RCP6),and one very high baseline emission scenario (RCP8.5).Global annual mean surface temperature at the end of the 21st century, compared to that at the end of the 20th century, was projected to rise by 0.3-1.7 • C under the RCP2.6 scenario and 2.6-4.8• C under the RCP8.5 scenario.At present, trends in greenhouse gas emissions are consistent with those assumed in the high emission scenario RCP8.5 [14].
The ecological impacts of global warming have been studied from polar terrestrial to tropical marine environments.It is evident that climate change can affect demography, geographic distribution, body size, fishery catch, phenology of populations, and biodiversity of ecosystems [15][16][17][18][19][20][21][22].Global meta-analyses indicated significant range shifts averaging 6.1 km per decade towards the poles (or meters per decade upward) and a significant mean advancement of spring events by 2.3 days per decade [23].Studies on the responses of ocean ecosystems to future climate change include projections on primary production [24], species distribution, and abundance at local [25,26] to global [27,28] scales under different emission scenarios.
Climate change in Japan due to global warming is inducing serious concerns among the Japanese public and political sectors [29,30].It was reported that the annual mean surface air temperature over Japan has increased at a rate of 1.19 • C/100 years, faster than the global average [31].Annual mean air temperatures in Japan were projected to rise by 1.1 ± 0.4 • C and 4.4 ± 0.6 • C under the RCP2.6 and RCP8.5 scenarios at the end of the 21st century, relative to the end of the 20th century [32].Although the warming impacts on the fishery species around Japan were mentioned in the report from the Ministry of Agriculture, Forestry and Fisheries [33], little information exists on the snow crab.Past warm phases provide a preview of the future, in which C. opilio in the Sea of Japan will shift northward and the resource decline will continue [34,35].
In this study, a survival-transport model was applied to clarify the impact of climate change on snow crab larvae in the Sea of Japan.We focus on future responses, but also summarize the changes simulated in the second half of the 20th century as a reference.We were not concerned with organismal-level physiological changes, individual-level behavioral changes, or ecosystem-level food-web-interaction changes, although all four were apparent and interlinked [36].
The paper is organized as follows: Section 2 introduces the model configuration.Section 3 describes changes in spawning and larval settlement in response to global warming.There is a brief discussion on the impacts of hydrodynamic factors in Section 4. Conclusions are presented at the end.

Model Descriptions
A survival and transport model was developed for snow crab larvae in the Sea of Japan (Figure 1a).The model included the process of releasing virtual particles from spawning grounds, the tracking and sinking processes during larval advection, accompanied by judging their survival status, and the final process of larval settlement on the seabed, with an additional diurnal vertical migration for the second zoea [6].The model can reproduce the distribution of crab larval drift and settlement on the seabed from 1999 to 2013 [37].Without regard for the impact of female crab abundance on the distribution of spawns, a uniform distribution of released particles was adopted as the initial condition in the simulation.Specifically, 25 particles were regularly distributed in the surface layer of the grid with a water depth ranging between 200 m and 500 m.If the SST of the grid fell within the range of 5 °C-16 °C, the particles were released as first zoea.The spawning sites in the Sea of Japan were divided into eight regions (Figure 1b).Regions A, B, and C were important bottom trawling fishery grounds for snow crab, respectively located west of Toyama Prefecture, west of Niigata Prefecture, and west of Hokkaido Prefecture [39].Region D was located east of the Korean peninsula.Region E was at Yamato Bank.Regions F and G were north of 46°N, along the east and west coast, respectively.The last spawning site, located between 40°N and 46°N along the west coast in the Sea of Japan, was named region H.
Previous studies have reported that snow crab spawning generally occurs in spring [6].However, we do not know how the hatching season would respond to global warming.It was previously reported that climate change accounted for 37% of a trend towards earlier egg-laying in breeding birds [40].Consequently, we assumed in our study that the spawning process could occur at any time of year, if SST was proper for first zoea spawning.To be specific, 37 release-drift events were simulated for a given year.The first event began on 1 January, with subsequent events initiated 10 days apart.Each event lasted for 120 days.The first zoea, second zoea, and megalopa stages were set at 30, 30, and 60 days, respectively [6,41,42].All hydrodynamic data were updated daily and were linearly interpolated into hourly frequencies using a time step of 3600 s.
During larval advection, the survival temperature limits were set at 5 °C-16 °C for the zoea stage (first 60 days in each cycle) and 5 °C-14 °C for the megalopa stage (second 60 days) [41].If one larva stayed out of the proper temperature range for 24 consecutive hours, it would "die" and was removed from the simulation.Two conditions were set for the larvae to settle on the seabed.The first required that the larvae had been transported for longer than 90 days, meaning that they have developed into the late megalopa stage and had prepared to settle on the seabed.The second The 3-D flow field was indispensable in the particle tracking process, while water temperature was necessary for the spawning process and survival judgement.In this study, we used the flow field and water temperature data from daily reanalysis datasets of Data Assimilation Research of the East Asian Marine System (DREAMS) [38] Without regard for the impact of female crab abundance on the distribution of spawns, a uniform distribution of released particles was adopted as the initial condition in the simulation.Specifically, 25 particles were regularly distributed in the surface layer of the grid with a water depth ranging between 200 m and 500 m.If the SST of the grid fell within the range of 5 • C-16 • C, the particles were released as first zoea.The spawning sites in the Sea of Japan were divided into eight regions (Figure 1b).Regions A, B, and C were important bottom trawling fishery grounds for snow crab, respectively located west of Toyama Prefecture, west of Niigata Prefecture, and west of Hokkaido Prefecture [39].Region D was located east of the Korean peninsula.Region E was at Yamato Bank.Regions F and G were north of 46 • N, along the east and west coast, respectively.The last spawning site, located between 40 • N and 46 • N along the west coast in the Sea of Japan, was named region H.
Previous studies have reported that snow crab spawning generally occurs in spring [6].However, we do not know how the hatching season would respond to global warming.It was previously reported that climate change accounted for 37% of a trend towards earlier egg-laying in breeding birds [40].Consequently, we assumed in our study that the spawning process could occur at any time of year, if SST was proper for first zoea spawning.To be specific, 37 release-drift events were simulated for a given year.The first event began on 1 January, with subsequent events initiated 10 days apart.Each event lasted for 120 days.The first zoea, second zoea, and megalopa stages were set at 30, 30, and 60 days, respectively [6,41,42].All hydrodynamic data were updated daily and were linearly interpolated into hourly frequencies using a time step of 3600 s.
During larval advection, the survival temperature limits were set at 5 • C-16 • C for the zoea stage (first 60 days in each cycle) and 5 • C-14 • C for the megalopa stage (second 60 days) [41].If one larva stayed out of the proper temperature range for 24 consecutive hours, it would "die" and was removed from the simulation.Two conditions were set for the larvae to settle on the seabed.The first required that the larvae had been transported for longer than 90 days, meaning that they have developed into the late megalopa stage and had prepared to settle on the seabed.The second condition required that the larvae be hyperbenthic and located at an area with a water depth between 200 m and 500 m.

Control Group
There were six numerical experiments of larval models in the control group, covering from past to future.They were all driven by the hydrodynamics of DREAMS datasets [38].DREAMS was simulated in three climate periods, the past, present, and the future.The future period consisted of data at the middle and end of the 21st century under the RCP2.6 and RCP8.5 scenarios.The air-sea boundary conditions for the 1960s and 2010s were from the JRA55 dataset [43,44].The change in air-sea boundary conditions for simulation in four future situations was a 1.0 • C increase of air temperature in RCP2.6 2050s, a 2.0 • C increase of air temperature in RCP8.5 2050s, a 1.0 • C increase of air temperature in RCP2.6 2100s, and a 4.0 • C increase of air temperature in RCP8.5 2100s.Therefore, there was a total of six subsets in the DREAMS datasets, each of which contained ten years of daily data (Table 1).The larval simulation in the control group was forced to the last five-year hydrodynamic data from each subset because of high uncertainty in the spin-up phase of DREAMS during the first several years.

Sensitive Groups
The sensitive simulation is structured to remove either the change in water temperature or that in flow field among different eras.Therefore, there are two sensitive groups, "Tmodern" and "UVmodern", each of which includes one past and four future experiments (Table 2).In the "Tmodern" group, all the experiments are forced by the same 2010s (referred to "modern") temperature (referred to as 'T') but a different flow field corresponding to one past and four future situations.On the other hand, all the experiments in group "UVmodern" are forced by different temperatures, but the same 2010s (referred to as "modern") flow field (referred to as "UV").As for the experiment name, "T" is short for water temperature; "U" and "V" are two horizontal components of velocity vector; and "mdn" is short for "modern" (2010s subset).Take "Tmdn RCP8.5 2050s" as example, it means the experiment was forced by water temperature from 2010s and flow field from the RCP8.5 2050s subset.

Climate Change Trend in DREAMS
The annual mean water temperature in the upper 200 m, where crab larvae reside, is presented in Figure A1.Water temperature in the Sea of Japan was found to consistently increase over time.It was estimated to warm by 0.6 • C from the 1960s to 2010s, and to continue rising by 1.4 • C and 2.8 • C under the RCP2.6 and RCP8.5 scenarios by the end of the 21st century.As for the flow field, a triple-branch pattern was observed in the flow field for the Tsushima Warm Current (TWC) (Figure A2a).The current path along Honshu Island represented the first branch of the TWC (FBTWC), the path along the continental shelf edge was the second branch (SBTWC), and the path along the east coast of the Korean peninsula represented the third branch (TBTWC) [45].SBTWC and TBTWC met at (132 • E, 37.5 • N), flowed northeastward, and then joined with FBTWC at (139 • E, 40 • N).The TWC in the past strengthened slightly for only two branches, FBTWC and SBTWC (Figure A2b).However, in the future group, the TWC strengthened by 30%-50%, especially for the TBTWC, which reached as far as 40 • N. As a result, the meeting location of SBTWC and TBTWC shifted to north of 38 • N and, afterwards, the joint current was also reinforced and shifted further north than in the present.Moreover, the northward flow off Hokkaido moved shoreward in the past and future when compared with the present.

Spawning of Snow Crab
Spawning numbers showed distinct seasonal variation, with much higher numbers from December to May than in other months and a peak in May (Figure 2b).In midsummer, spawning can only occur in regions F, G, and H, at quite low levels (Figure 2f,h).Two types of seasonal variation could be distinguished in the spawning sites.The variations in spawning sites south of 41. 5 • N (regions A, B, D,  and E) demonstrated that particles were released at the maximum level from January to May and then dramatically decreased until November (Figure 2, left panel).However, at the sites north of 41.5 • N (regions C, F, G, and H), two spawning peaks were identified in May-June and October-November (Figure 2, right panel).
The seasonal variation in crab spawning was fully dependent on the spatio-temporal variations of SST.We therefore classified the possible spawning sites identified in Figure 1b by SST as always/sometimes/never falling in the range of 5 • C-16 • C (Figure 3).In the 2010s, the highest spawning levels occurred in May, as indicated in Figure 3 by the limited blue area that month.In August, the particles could only be released in regions G and H, where SST sometimes fell in the proper range.The two spawning peaks in regions north of 41.5 • N could be distinguished easily due to the larger area covered in red during May and November.The seasonal variation in crab spawning was fully dependent on the spatio-temporal variations of SST.We therefore classified the possible spawning sites identified in Figure 1b by SST as always/sometimes/never falling in the range of 5 °C-16 °C (Figure 3).In the 2010s, the highest spawning levels occurred in May, as indicated in Figure 3 by the limited blue area that month.In August, the particles could only be released in regions G and H, where SST sometimes fell in the proper range.The two spawning peaks in regions north of 41.5°N could be distinguished easily due to the larger area covered in red during May and November.When comparing seasonal variations among experiments, minimal decline was found in the regions south of 41.5°N in the future (Figure 2, left panel).In contrast, spawning in region C increased remarkably in future spring, with larger increments under the RCP8.5 scenario (Figure 2d).In regions F, G, and H (Figure 2f and h), the spawning peak in spring shifted forward, while that in autumn shifted backward, showing the phenological effects of warming [19].For example, the When comparing seasonal variations among experiments, minimal decline was found in the regions south of 41.5 • N in the future (Figure 2, left panel).In contrast, spawning in region C increased remarkably in future spring, with larger increments under the RCP8.5 scenario (Figure 2d).In regions F, G, and H (Figure 2f,h), the spawning peak in spring shifted forward, while that in autumn shifted backward, showing the phenological effects of warming [19].For example, the spring peak in region H occurred on 10 June in the 1960s and was advanced 10 days in the 2010s and another 10 days/30 days in the 2100s, under the RCP2.6/RCP8.5 scenarios (Figure 2h).
Figure 4 presents the quantitative estimations of spawning amounts in the six experiments.As the ocean becomes warmer in the future (Figure A1), the area of SST > 16 • C in the southern region would be larger, as would the area of SST > 5 • C in the northern region.This resulted in different spawning trends between regions south and north of 41.5 • N. Additionally, the variation in spawning amounts was more significant under the RCP8.5 scenario than RCP2.6.For example, the spawning amount in region A was reduced by 23% in RCP8.5 2100s, but 6% in RCP2.6 2100s.When comparing seasonal variations among experiments, minimal decline was found in the regions south of 41.5°N in the future (Figure 2, left panel).In contrast, spawning in region C increased remarkably in future spring, with larger increments under the RCP8.5 scenario (Figure 2d).In regions F, G, and H (Figure 2f and h), the spawning peak in spring shifted forward, while that in autumn shifted backward, showing the phenological effects of warming [19].For example, the spring peak in region H occurred on 10 June in the 1960s and was advanced 10 days in the 2010s and another 10 days/30 days in the 2100s, under the RCP2.6/RCP8.5 scenarios (Figure 2h).
Figure 4 presents the quantitative estimations of spawning amounts in the six experiments.As the ocean becomes warmer in the future (Figure A1), the area of SST > 16 °C in the southern region would be larger, as would the area of SST > 5 °C in the northern region.This resulted in different spawning trends between regions south and north of 41.5°N.Additionally, the variation in spawning amounts was more significant under the RCP8.5 scenario than RCP2.6.For example, the spawning amount in region A was reduced by 23% in RCP8.5 2100s, but 6% in RCP2.6 2100s.Unexpectedly, the spawning amount in the 1960s was higher than today in most regions.The seasonal variation depicted in Figure 2   Unexpectedly, the spawning amount in the 1960s was higher than today in most regions.The seasonal variation depicted in Figure 2

Larval Settlement
After release, the crab larvae were advected by the cyclone-like circulation in the Sea of Japan, mainly consisting of the Liman Current, the North Korea Cold Current, the Ulleung Eddy, and TWC

Larval Settlement
After release, the crab larvae were advected by the cyclone-like circulation in the Sea of Japan, mainly consisting of the Liman Current, the North Korea Cold Current, the Ulleung Eddy, and TWC [46], and sank.In the later megalopa stage, some larvae settled on the seabed off Honshu and Hokkaido Islands (regions A, B, and C; Figure 6).The fewest number of larvae settled in the 2010s.The settling rate, i.e., the ratio of settled larvae number to initially released number, was also the lowest one in the 2010s.Among the settlement, 99% occurred in regions A, B, and C.However, in RCP8.5 2100s, this percent decreased to 88.5% because some megalopa settled on the seabed north of 46 • N in region F, while others settled along the coast in region H and a few settled in region E (Figure 6d).Additionally, the gathering site of settled megalopa shifted northward as climate changed.The highest settlement proportion occurred in region A in the 1960s, but shifted to region B in the 2010s.From then on, region B exhibited a steady decline.As a result, the settlement proportion in region C increased by 27%(RCP2.6)-45%(RCP8.5) in the 2100s.The effect of different RCP scenarios on crab larval settlement was much more evident in Figure 7. Settlement in regions A and B was lower under the RCP8.5 scenario than RCP2.6, in response to much warmer ocean temperatures under the high GHG emission scenario.However, settlement was highest in RCP8.5 2100s, mainly due to settlement in region C. Region C was also responsible for similarly high levels of total settlement in the 2050s, under two different RCPs.However, settlement in region F and other regions was only evident in RCP8.5 2100s.The effect of different RCP scenarios on crab larval settlement was much more evident in Figure 7. Settlement in regions A and B was lower under the RCP8.5 scenario than RCP2.6, in response to much warmer ocean temperatures under the high GHG emission scenario.However, settlement was highest in RCP8.5 2100s, mainly due to settlement in region C. Region C was also responsible for similarly high levels of total settlement in the 2050s, under two different RCPs.However, settlement in region F and other regions was only evident in RCP8.5 2100s.In addition to the departure from the Sea of Japan through the Tsugaru Strait and the Soya Strait, four final states of the larvae were identified, i.e., death due to low or high temperatures, drifting in the ocean, or settlement on the seabed (Table 3).In the future experiments, increased settlement resulted from a decline in the low-temperature (< 5 °C) mortality rate.A 1%-3% decrease in mortality can result in a 1-2 time increase in settlement amount.Spatial distribution analysis indicated that low-temperature mortality drastically decreased in the southern Sea of Japan (typical area: 134-136°E, 36-38°N) and in northwest Hokkaido in the future experiments (Figure 8).For example, as shown in Figure A3, water with temperature > 5 °C in the southern Sea of Japan was more abundant in RCP2.6 2050s than in the 2010s.A similar discrepancy in water temperature was observed for low-temperature mortality distribution (Figure A3c and Figure 8e).In addition to the departure from the Sea of Japan through the Tsugaru Strait and the Soya Strait, four final states of the larvae were identified, i.e., death due to low or high temperatures, drifting in the ocean, or settlement on the seabed (Table 3).In the future experiments, increased settlement resulted from a decline in the low-temperature (<5 • C) mortality rate.A 1%-3% decrease in mortality can result in a 1-2 time increase in settlement amount.Spatial distribution analysis indicated that low-temperature mortality drastically decreased in the southern Sea of Japan (typical area: 134-136 • E, 36-38 • N) and in northwest Hokkaido in the future experiments (Figure 8).For example, as shown in Figure A3, water with temperature >5 • C in the southern Sea of Japan was more abundant in RCP2.6 2050s than in the 2010s.A similar discrepancy in water temperature was observed for low-temperature mortality distribution (Figures A3c and 8e).

Hydrodynamic Factors Controlling Settlement
As mentioned in Section 2.2.2, the "Tmodern" and "UVmodern" groups were run to separate the effects of changes in water temperature or changes in flow field.The difference in settlement between "Tmodern" and "control" was more distinct than that between "UVmodern" and "control" (Figure 9), suggesting the dominant impact of water temperature in the Sea of Japan.This influence was observed in all periods from the past to the future.For example, the annual settlement in region B, under scenario of RCP2.6 2100s, was reduced from 1740 in "control" to 322 in "Tmodern", but increased to 2093 in the "UVmodern" group.This indicated that a 1.4 °C increase in water temperature under the high mitigation scenario of RCP2.6 would increase larval settlement by five times and reveals the rectification of flow field on the distribution of larval settlement.As Figure A2

Hydrodynamic Factors Controlling Settlement
As mentioned in Section 2.2.2, the "Tmodern" and "UVmodern" groups were run to separate the effects of changes in water temperature or changes in flow field.The difference in settlement between "Tmodern" and "control" was more distinct than that between "UVmodern" and "control" (Figure 9), suggesting the dominant impact of water temperature in the Sea of Japan.This influence was observed in all periods from the past to the future.For example, the annual settlement in region B, under scenario of RCP2.6 2100s, was reduced from 1740 in "control" to 322 in "Tmodern", but increased to 2093 in the "UVmodern" group.This indicated that a 1.4 • C increase in water temperature under the high mitigation scenario of RCP2.6 would increase larval settlement by five times and reveals the rectification of flow field on the distribution of larval settlement.As Figure A2 shows, FBTWC in the 2010s was stronger than that in the RCP2.6 2100s, with a southerly weaker joint current of SBTWC and TBTWC.Therefore, the settlement in region B was a little higher in "UVmodern" than in "control".shows, FBTWC in the 2010s was stronger than that in the RCP2.6 2100s, with a southerly weaker joint current of SBTWC and TBTWC.Therefore, the settlement in region B was a little higher in "UVmodern" than in "control".To elucidate the effect of flow field on settlement, we identified the settlement proportion of five regions in four future cases from the "control" and "UVmodern" groups (Table 4).When using the flow field of the 2010s, settlement in region B was higher, as mentioned above, while settlement in the region north of 41.5°N (regions C and F) was lower.This was related to the offshore shift and weaker flow off west Hokkaido in the 2010s (Figure A2).Although the variation due to the change in flow field was distinguishable, it was ~14% at most, much weaker than the impact of water temperature.However, RCP8.5 2100s differed from the others, as settlement in the region north of 41.5°N was 3% higher in "UVmodern" than in "control".To elucidate the effect of flow field on settlement, we identified the settlement proportion of five regions in four future cases from the "control" and "UVmodern" groups (Table 4).When using the flow field of the 2010s, settlement in region B was higher, as mentioned above, while settlement in the region north of 41.5 • N (regions C and F) was lower.This was related to the offshore shift and weaker flow off west Hokkaido in the 2010s (Figure A2).Although the variation due to the change in flow field was distinguishable, it was ~14% at most, much weaker than the impact of water temperature.However, RCP8.5 2100s differed from the others, as settlement in the region north of 41.5 • N was 3% higher in "UVmodern" than in "control".

Clarification of the Simulation
The results of six experiments in the control group suggested a positive impact of global warming on crab larval settlement in the Sea of Japan.These results contrasted with the declining trends predicted for Newfoundland and Labrador snow crab under a warming oceanographic regime [35].However, the projected increase in larval settlement does not necessarily signify an increase in snow crab resources in the future.Predation by jellyfish and juvenile pollock on snow crab planktonic larvae was not included in the present model.This may increase under warming temperatures and diminish larval settlement in the future.Though the warming conditions had a positive effect on crab size and a negative effect on productivity and recruitment [47], a strong negative relationship was present between the bottom temperature and early-life survival of small crabs [48].As indicated in Figure A4, the bottom temperature in the Sea of Japan was projected to rise by 3 • C at the end of the 21st century, under the RCP8.5 scenario, which would present an extreme challenge for young benthic crabs.Moreover, predation of small crabs by benthic fish, as well as cannibalism by elder crabs, can lead to the mortality of metamorphosed benthic crabs.Since it takes 7-8 years for crabs to develop to maturity, the effects of climate change on the crab abundance are unclear, although the harvesting of snow crab off Hokkaido Island is likely to increase in the future.

Conclusions
In this paper, the change in crab spawning and larval settlement as the climate changes is analyzed with help of the survival and transport model.Seasonal variation in spawning was distinct in the Sea of Japan and closely related to SST.Spawning in regions south of 41.5 • N was maximal from January to May, and then dramatically decreased until November.The annual spawning level at these sites is predicted to decline slightly in the future.In the regions north of 41.5 • N, two spawning peaks occur every year.In response to climate change, the spring peak would be advanced and the autumn peak delayed.The annual spawning amount would increase in the future, with larger increments under the RCP8.5 scenario.
We found that ~90% of larval settlement occurred in regions A, B, and C, the most important benthic crab grounds along the Japanese coast.Settling in the 2010s was the lowest among the six experiments.The reduction of low-temperature mortality in crab larvae was found to be the primary reason for settlement increase in the future.Additionally, the gathering site of settled megalopa likely shifted northward due to global warming, with greater contribution from settling in region C. We further investigated the impacts of change in either water temperature or flow fields on the larval settlement.Water temperature was found to play a predominant role in larval settlement, not only in the amount, but also in the distribution.However, the difference in flow field between the future and the present could modulate the settlement off Hokkaido.
In this study, the change in snow crab larvae in response to global warming was reflected by the balance among rates of mortality, drifting, and settlement on the seabed.Therefore, this study linked changes in snow crab larvae to changes at the population level, the third interlinked level of biological organization [36].The model has several shortcomings.One limitation is that the duration of the zoea I, zoea II, and megalopa stages was fixed in the model.The duration of these stages may change due to the warming ocean.Another weakness is that consumption and predation of crab larva were not considered.Plankton blooms are thought to be accelerated due to the phenological effects of global warming [23] and, thereby, cause changes in the consumption and growth of crab larvae.To improve the crab larval model, a lower trophic level ecosystem model should be coupled, which will lead to a more accurate simulation.Finally, it needs to be reiterated that the projected increase in larval settlement does not signify an increase of snow crab resources in the future, for many factors can lead to the mortality of metamorphosed benthic crabs.However, we believe that the snow crab fishery off Hokkaido Island will be more significant than today.
Funding: This research was funded by the Environment Research and Technology Development Fund of the Ministry of the Environment, Japan, grant number S-13; the National Natural Science Foundation of China, grant number 41576010.larval settlement.Water temperature was found to play a predominant role in larval settlement, not only in the amount, but also in the distribution.However, the difference in flow field between the future and the present could modulate the settlement off Hokkaido.
In this study, the change in snow crab larvae in response to global warming was reflected by the balance among rates of mortality, drifting, and settlement on the seabed.Therefore, this study linked changes in snow crab larvae to changes at the population level, the third interlinked level of biological organization [36].The model has several shortcomings.One limitation is that the duration of the zoea I, zoea II, and megalopa stages was fixed in the model.The duration of these stages may change due to the warming ocean.Another weakness is that consumption and predation of crab larva were not considered.Plankton blooms are thought to be accelerated due to the phenological effects of global warming [23] and, thereby, cause changes in the consumption and growth of crab larvae.To improve the crab larval model, a lower trophic level ecosystem model should be coupled, which will lead to a more accurate simulation.Finally, it needs to be reiterated that the projected increase in larval settlement does not signify an increase of snow crab resources in the future, for many factors can lead to the mortality of metamorphosed benthic crabs.However, we believe that the snow crab fishery off Hokkaido Island will be more significant than today.
Author Contributions: Conceptualization, Mao X. and Guo X.; methodology, Wang Y.; validation, Wang Y. and X.; writing, Mao X. and Guo X.; calculation of DREAMS, Takayama K.

Sustainability 2019 , 19 Figure 1 .
Figure 1.Model domain.(a) Bathymetry in the Sea of Japan; (b) Possible spawning sites of snow crab, determined by water depth in the range of 200 m to 500 m.

Figure 1 .
Figure 1.Model domain.(a) Bathymetry in the Sea of Japan; (b) Possible spawning sites of snow crab, determined by water depth in the range of 200 m to 500 m.
. The DREAMS system has two sub-models with a low resolution, one covering almost the entire Pacific Ocean and a high resolution one covering the northwestern Pacific Ocean.The low-resolution sub-model was run for 150 years from 1960 to 2100, whose results provide the lateral boundary conditions and initial conditions for the high-resolution one.The horizontal resolution of the high-resolution sub-model was 5' in longitude and 4' in latitude.There were 40 layers in the vertical direction (1 m, 3.5 m, 6.5 m, 10 m, 15 m, 22 m, 30 m, 39 m, 50 m, 64 m, 81 m, 100 m, 120 m, . . ., 5674 m).

Figure 2 .
Figure 2. Annual time-series of spawning levels in six experiments of the control group.

Figure 2 . 19 Figure 3 .
Figure 2. Annual time-series of spawning levels in six experiments of the control group.Sustainability 2019, 11, x FOR PEER REVIEW 7 of 19

Figure 3 .
Figure 3. Seasonal variation of the SST range classification in 2010s.(a) February; (b) May; (c) August; and (d) November.Red dots indicate that the SST in those grids always fell within the range of 5 • C-16 • C; brown indicates the SST may have been out of range in some years, and blue means the SST never fell in the range during 2009-2013.

Figure 3 .
Figure 3. Seasonal variation of the SST range classification in 2010s.(a) February; (b) May; (c) August; and (d) November.Red dots indicate that the SST in those grids always fell within the range of 5 °C-16°C; brown indicates the SST may have been out of range in some years, and blue means the SST never fell in the range during 2009-2013.

Figure 4 .
Figure 4. Annual spawning amount in different regions among six experiments.
indicated a distinct increase in November in regions A and B and in August-September in regions F and G.A comparison of SST range classifications between the 2010s and 1960s (Figure 3c,d and Figure 5a,b) provides clarification for these unexpected results.More particles were released from region F in August and from region A in November in the 1960s than in the present period.Sustainability 2019, 11, x FOR PEER REVIEW 8 of 19

Figure 4 .
Figure 4. Annual spawning amount in different regions among six experiments.
indicated a distinct increase in November in regions A and B and in August-September in regions F and G.A comparison of SST range classifications between the 2010s and 1960s (Figure 3c-d and Figure 5a-b) provides clarification for these unexpected results.More particles were released from region F in August and from region A in November in the 1960s than in the present period.

Figure 5 .
Figure 5. Same as Figure 3, but for the 1960s.(a) August and (b) November.

Figure 5 .
Figure 5. Same as Figure 3, but for the 1960s.(a) August and (b) November.

Figure 6 .
Figure 6.Distribution of total settled larvae (shading) in six experiments and zonal statistics.The percent after the settled number in the top left is the settling rate; the percent after the region is the settlement proportion of that region.(a) 2010s; (b) 1960s; (c) RCP8.5 2050s; (d) RCP8.5 2100s; (e) RCP2.6 2050s; and (f) RCP2.6 2100s.

Figure 6 .
Figure 6.Distribution of total settled larvae (shading) in six experiments and zonal statistics.The percent after the settled number in the top left is the settling rate; the percent after the region is the settlement proportion of that region.(a) 2010s; (b) 1960s; (c) RCP8.5 2050s; (d) RCP8.5 2100s; (e) RCP2.6 2050s; and (f) RCP2.6 2100s.

Figure 7 .
Figure 7. Annual settlement amount in different regions among six experiments.

Figure 7 .
Figure 7. Annual settlement amount in different regions among six experiments.

Figure 8 .
Figure 8. Distribution of low-temperature mortality (shading) in the 2010s (a) and the difference between other experiments and the 2010s (b-f).Percent in (a) is the death rate from cold in 2010s; percent in (b-f) shows the change in death rate from cold from the 2010s to the other experiments.

Figure 8 .
Figure 8. Distribution of low-temperature mortality (shading) in the 2010s (a) and the difference between other experiments and the 2010s (b-f).Percent in (a) is the death rate from cold in 2010s; percent in (b-f) shows the change in death rate from cold from the 2010s to the other experiments.

Figure 9 .
Figure 9.Comparison of annual settlement amounts among three groups.(a) Control group; (b) Tmodern group; and (c) UVmodern group.

Figure 9 .
Figure 9.Comparison of annual settlement amounts among three groups.(a) Control group; (b) Tmodern group; and (c) UVmodern group.

Figure 2 .
Figure 2. Climatological monthly mean flow field at 50 m in August in six subsets.The green line represents the isobath of 500 m.The color in (a) denotes the flow magnitude (cm/s) in the 2010s; the color in (b-f) denotes the difference in flow magnitude between the other subsets and the 2010s.

Figure A2 .
Figure A2.Climatological monthly mean flow field at 50 m in August in six subsets.The green line represents the isobath of 500 m.The color in (a) denotes the flow magnitude (cm/s) in the 2010s; the color in (b-f) denotes the difference in flow magnitude between the other subsets and the 2010s.Sustainability 2019, 11, x FOR PEER REVIEW 16 of 19

Figure 3 .
Figure 3. Climatological monthly mean SST in August in the 2010s RCP2.6 2050s (b), and the difference (c).The green line in (c) shows the isobath of 500 m.

Figure A3 .
Figure A3.Climatological monthly mean SST in August in the 2010s (a), RCP2.6 2050s (b), and the difference (c).The green line in (c) shows the isobath of 500 m.

Figure 3 .
Figure 3. Climatological monthly mean SST in August in the 2010s (a), RCP2.6 2050s (b), and the difference (c).The green line in (c) shows the isobath of 500 m.

Figure A4 .
Figure A4.Climatological monthly mean temperature at the bottom layer (water depth <1000 m) in June in six subsets.(a) 2010s; (b-f) the differences between the other subsets and the 2010s.

Table 1 .
Name of larval simulation in the control group and relevant DREAMS information.

Table 2 .
Experiments in two sensitive groups and relevant subset of hydrodynamic data.

Table 3 .
The number of released (initial) and ending particles in six experiments.

Table 3 .
The number of released (initial) and ending particles in six experiments.

Table 4 .
Settlement proportion in different regions between "control" and "UVmodern".Results of control group are shown by grey background.