Continuous Monitoring and Future Projection of Ocean Warming, Acidification, and Deoxygenation on the Subarctic Coast of Hokkaido, Japan

As the ocean absorbs excessive anthropogenic CO2 and ocean acidification proceeds, it is thought to be harder for marine calcifying organisms, such as shellfish, to form their skeletons and shells made of calcium carbonate. Recent studies have suggested that various marine organisms, both calcifiers and non-calcifiers, will be affected adversely by ocean warming and deoxygenation. However, regardless of their effects on calcifiers, the spatiotemporal variability of parameters affecting ocean acidification and deoxygenation has not been elucidated in the subarctic coasts of Japan. This study conducted the first continuous monitoring and future projection of physical and biogeochemical parameters of the subarctic coast of Hokkaido, Japan. Our results show that the seasonal change in biogeochemical parameters, with higher pH and dissolved oxygen (DO) concentration in winter than in summer, was primarily regulated by water temperature. The daily fluctuations, which were higher in the daytime than at night, were mainly affected by daytime photosynthesis by primary producers and respiration by marine organisms at night. Our projected results suggest that, without ambitious commitment to reducing CO2 and other greenhouse gas emissions, such as by following the Paris Agreement, the impact of ocean warming and acidification on calcifiers along subarctic coasts will become serious, exceeding the critical level of high temperature for 3 months in summer and being close to the critical level of low saturation state of calcium carbonate for 2 months in mid-winter, respectively, by the end of this century. The impact of deoxygenation might often be prominent assuming that the daily fluctuation in DO concentration in the future is similar to that at present. The results also suggest the importance of adaptation strategies by local coastal industries, especially fisheries, such as modifying aquaculture styles.


INTRODUCTION
Global change-driven global ocean warming, acidification, and deoxygenation are all thought to be caused primarily by excessive anthropogenic CO 2 emissions. Various studies have suggested that the rise in water temperature caused by global change will alter habitats and marine ecosystem functions [e.g., Intergovernmental Panel on Climate Change (IPCC), 2014]. However, the impacts of ocean acidification and deoxygenation are difficult to detect directly, and relatively little evidence has been reported.
Ocean acidification leads to lower pH in the surface ocean and the pH has already been lowered by 0.1 globally over the past century (Orr et al., 2005). The global average oceanic pH is currently around 8.1 and is anticipated to decrease to 7.8 by the end of this century (e.g., Feely et al., 2008; Intergovernmental Panel on Climate Change (IPCC), 2014). Recent studies have suggested that ocean acidification will affect marine organisms that form shells or bodies made of calcium carbonate (calcifiers) because these chemicals are more easily dissolved or harder to generate in high-CO 2 , low-pH environments. The impact of ocean acidification on calcifiers has already been reported, such as the significant mortality of oyster (Crassostrea gigas) larvae (Feely et al., 2008;Barton et al., 2012) and damage to the carapaces of larval Dungeness crabs (Metacarcinus magister; Bednaršek et al., 2020) in the high-CO 2 environment of the Pacific northwest coast of the United States of America. Several studies have shown that ocean acidification also affects the life stages of calcifiers differently, and the larvae are more vulnerable to high CO 2 (e.g., Kurihara, 2008;Kimura et al., 2011;Waldbusser et al., 2015;Zhan et al., 2016;Onitsuka et al., 2018;Foo et al., 2020). Ocean acidification also affects non-calcifiers, such as by impairing the fertilization of eggs (e.g., Morita et al., 2009) and olfactory discrimination and homing ability (e.g., Munday et al., 2009Munday et al., , 2014Dixson et al., 2015). Therefore, similar to ocean warming, ocean acidification also affects marine ecosystems widely.
Deoxygenation has long been recognized as a local phenomenon, such as that caused by blooms of red tides and subsequent extremely low dissolved oxygen (DO) concentration environments in the bottom layers. Recently, world-wide deoxygenation has been observed as a result of ocean warming followed by lower oxygen saturation in warmer waters, reduced air-sea oxygen exchange with stratification of the surface layer and increased consumption of DO by enhanced respiration and dissolution of organic matter at higher temperatures (e.g., Turner et al., 2005;Vaquer-Sunyer and Duarte, 2008;Rabalais et al., 2010;Frieder et al., 2012;Marumo and Yokota, 2012;Schmidtko et al., 2017;Christian and Ono, 2019). Deoxygenation affects various marine organisms more directly and acutely than ocean acidification. Moreover, recent studies have suggested combined multi-stressor synergistic effects of ocean acidification and deoxygenation on calcifiers that may magnify the adverse impacts (e.g., Melzner et al., 2013;DePasquale et al., 2015;Gobler and Baumann, 2016;IPCC, 2018).
Recent observational studies have shown that both physical and biogeochemical parameters regulating ocean acidification and deoxygenation vary markedly in space and time. For example, the spatiotemporal variability is generally much larger in coastal oceans than in open oceans (e.g., Wootton et al., 2008;Dai et al., 2009;Rabalais et al., 2010;Hofmann et al., 2011;Frieder et al., 2012;Baumann et al., 2014). The greater temporal variability is characterized by the larger seasonal and daily variation in coastal ocean parameters. Daily minimum saturation state value of aragonite (a polymorph of calcium carbonate; arag ) of around 1, which is close to the undersaturated and is unsuitable for calcifiers, have already been recorded at night in some coastal oceans (e.g., Frieder et al., 2012). The daily minimum DO is similarly low in the coastal oceans, where deoxygenation is caused by coexisting drivers of global change and local human impacts, such as additive oxygen utilization stimulated by nutrient inputs from the land and aquaculture (e.g., Melzner et al., 2013;Wallace et al., 2014).
However, fewer studies have examined the effects of ocean acidification and deoxygenation on subarctic coasts compared to tropical and subtropical coasts (e.g., Hofmann et al., 2011), although several have examined the subarctic open ocean (e.g., Wakita et al., 2017). This is mainly because of the difficulty in deploying multi-year monitoring, especially in winter when bad weather often prevents observations. It is important to elucidate the impact of global change, especially ocean acidification, on subarctic coasts for two reasons. First, ocean acidification is considered to affect calcifiers in subarctic oceans more significantly than in tropical, subtropical, and temperate oceans, because high-latitude systems are generally less buffered than low latitude systems and likely experience a higher frequency occurrence of under-saturated conditions. The other is because there is a high local fishery dependence on marine calcifiers. Most of the subarctic coasts in Japan are in Hokkaido Prefecture, which also accounts for one-fourth of the total fish catch, half of the calcifier catch, and 60% of the shellfish catch in Japan (Ministry of Agriculture, Forestry and Fisheries website: https://www.maff. go.jp/e/index.html). Therefore, it is very important to elucidate the current status of ocean acidification and deoxygenation and predict the possible future status from the perspective of local societies and the marine ecosystem of the subarctic coasts.
To assess the impacts of ocean acidification and deoxygenation on coastal organisms properly, we need to make long-term, high-frequency measurements of environmental parameters in subarctic regions. To investigate the diurnal variation in physical and biogeochemical properties related to ocean warming and acidification and deoxygenation in the subarctic region, this study measured physical and biogeochemical parameters to elucidate the temporal variation and its causes in the subarctic coastal region. The results were used for future projection toward the end of this century. Based on the combined results, we hope to develop scientific guidelines to combat the impact of global change on marine ecosystems and relevant local industries along subarctic coasts.

Monitoring Site
We established a site for monitoring and modeling physical and biogeochemical properties in Oshoro Bay, Hokkaido, Japan, as a representative subarctic coast site in Japan (Figure 1; Christian and Ono, 2019; Yamaka, 2019). The site is in the subarctic region of the Japan Sea and the longitude and latitude is 140.86 • E and 43.21 • N, respectively. Fisheries are a main local industry, and the target species include sea urchin (Strongylocentrotus intermedius), kelps (Laminaria spp.), Ezo abalone (Haliotis discus hannai), and scallops (Mizuhopecten yessoensis); all are representative subarctic fisheries species (Yamaka, 2019).

Monitoring
Since March 30, 2017, we have measured continuously temperature and salinity as physical parameters, and DO and pH (total scale) as biogeochemical parameters. A conductivity and temperature sensor (INFINITY-CTW ACTW-USB; JFE Advantech) was used to measure temperature and salinity hourly, while DO was measured hourly with a RINKO W AROW-USB (JFE Advantech). The accuracy of temperature, conductivity and DO is ±0.01 • C, ±0.01 mS cm −1 , and non-linearity ± 2% FS (at 1 atm, 25 • C), respectively. The range of the DO sensor was 0-20 mg l −1 DO concentration (or 0-200% DO saturation) and data beyond this range were eliminated. Calibration of the DO sensor was carried out by two-point (zero and span) calibration using 0 and 100% (saturated) oxygen waters. The two sensors have mechanical wipers that periodically sweep the sensing surface to keep the initial accuracy, avoiding bio-fouling without chemical materials and frequent maintenances.
To measure pH hourly, glass electrode pH sensors (SP-11 and SPS-14; Kimoto Electric) were used. The reproducibility of pH sensors is ±0.001. The sensors were calibrated against two seawater scale pH standard buffer solutions, 2-amino-2hydroxymethyl-1, 3-propanediol (TRIS) and 2-aminopyridine (AMP). The sensors were withdrawn for maintenance and calibration at 1-to 3-month intervals to minimize measurement biases caused by bio-fouling and measurement drift. At each calibration, it was confirmed that the sensor kept its Nernst response (ca. 0.05916V/pH at 25 • C and 35.0 salinity) within the range of ±0.5%.
Seawater samples were collected when the sensors were calibrated and analyzed using a total alkalinity (TA) titration analyzer (ATT-05; Kimoto Electric) to obtain the TA and a coulometer (MODEL 3000A; Nippon ANS) to determine the dissolved inorganic carbon (DIC) concentration (Wakita et al., 2017(Wakita et al., , 2021Yamaka, 2019). The values were calibrated against certified reference material (CRM) provided by Prof. A. G. Dickson (Scripps Institution of Oceanography, University of California San Diego) and KANSO. The precision of the analyses of both the DIC and TA was smaller than ±2 µmol kg −1 based on replicate samples. Then, the pH (total scale) at the in situ temperature were calculated from the carbonate dissociation constants in Lueker et al. (2000) and temperature, salinity, TA, and DIC using CO2SYS (Pierrot et al., 2006). pH at time t [pH(t)] obtained by the sensor was then corrected its drift by the following equation: where pH m (t) represents measured value of pH by the sensor at the time t, pH sample (t e ) and pH m (t e ) is the pH at the end time of each deployment (t e ) obtained by the seawater sample and sensor, respectively, and t i is the start time of each deployment.
Equation (1) assumes that drift of the pH sensor increases linearly with time over the whole period of single deployment. However, actual time course of sensor drift is not linear, so this may cause significant error in the pH(t) value after the drift correction. In 2020, we made a 55-days deployment of the same FIGURE 2 | Time-series of pH values in the Miyako Bay in August through October 2020. Dots in blue denote raw pH values obtained by a pH sensor (SP-11; pH raw ). Those in red demote pH values calculated by referring to pH values of water samples (pH sample ). Those in light blue and green show pH values after calibration and bias correction by pH values obtained by water samples weekly (pH 1 ; in light blue) and bi-monthly (at the initial and end time of the deployment; pH 2 ; in green).
pH sensor (SP-11) that we deployed in this study into the Miyako Bay, a bay in northern Japan having similar feature with that of the Oshoro Bay. At this time, we calculated pH(t) using the same Eq. (1) but with pH sample data obtained by weekly basis [pH 1 (t)] and those obtained by bi-monthly basis [at the initial and end time of the deployment; pH 2 (t)], and compared the drift-corrected values (Figure 2). We found that the root mean square of difference between pH 1 (t) and pH 2 (t) over the whole deployment period was 0.015. We interpreted this value as the measure of the uncertainty of pH(t) originated from different frequency of the drift correction (10 days and 2 months in this case). Therefore, we evaluate that the uncertainty of pH(t) calculated in this study is in the same order and is small enough for monitoring in the coastal regions with relatively larger daily and seasonal variations. arag of the seawater samples were calculated from the carbonate dissociation constants in Lueker et al. (2000) and temperature, salinity, TA, and DIC using CO2SYS (Pierrot et al., 2006). arag was also estimated using CO2SYS, with TA, the monitored temperature, and salinity. As no continuous monitoring TA data exist, we used the TA value obtained with a regression equation for north of 30 • N in the North Pacific (Lee et al., 2006): where the longitude of the study site is 140.86 • E ( Figure 1B).

Modeling
The current and future states of physical and biogeochemical parameters were simulated using the Regional Ocean Modeling System (ROMS). Of the versions of ROMS, we chose ROMS-Agrif (Penven et al., 2006), which embeds an ecosystem model (Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES); Aumont et al., 2003) in the ROMS. The bathymetry was expressed using the Earth 1 Arc Minute Topography Model (ETOPO1; Amante and Eakins, 2009). The ROMS was downscaled (nested) twice ( Figure 1A; Yamaka, 2019;Akamatsu, 2020). The first ROMS (L1) had a horizontal resolution of 1/3 degree and a temporal resolution of 1 h. To simulate physical and biogeochemical parameters using L1, boundary and initial conditions were needed. For this simulation, the Comprehensive Ocean-Atmosphere Data Set (COADS) 2005 (Da Silva et al., 1994) was used for the atmospheric forcing. For the oceanic physical and biogeochemical forcing, the World Ocean Atlas (WOA) 2009 (Goyet et al., 2000;Aumont and Bopp, 2006;Antonov et al., 2010;Garcia et al., 2010;Locarnini et al., 2010) was used. For future projection, the climate model outputs of the MIROC-EMS (model description and basic results of CMIP5-20c3m; Watanabe et al., 2011). Representative Concentration Pathways (RCP) 2.6 and 8.5 scenarios were used. The outputs of L1 were used as boundary and initial conditions to drive a second ROMS (L2) with a horizontal resolution of 1/15 degree and a temporal resolution of 1 h. The model was spun up for 5 years by driving with the 10-year mean current (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) and future (2086-2095) boundary conditions and the model outputs for the last 5 years were averaged and analyzed in this study.
When the monthly-averaged model results were compared with the corresponding observed results, the model results still differed from the observed results and the extent of modelobservation misfits differs largely among parameters and months (Figure 3). One of the possible reasons for the model-observation misfits is because the spatial resolution of L2 is still a bit too coarse to simulate the Oshoro bay. The other possible reason is because the model does not embed local processes such as inputs of freshwater and terrestrial nutrients to the coast, and primary production by seagrass/seaweeds. For example, as explained in 3.1, substantially low salinity was observed from late winter to early spring (in March and April), which was caused by the inflow of snow-melt water discharged directly to the bay from the surrounding land (Nakata et al., 2001). Such observed processes were not able to be reproduced by the model. Also, as well as phytoplankton bloom, abrupt primary production by seagrass/seaweeds is dominant in spring (in April through June) in the study site, which is not reproduced explicitly by the ecosystem model used in this study in which phytoplankton are assumed to be primary producers. This might lead to relative underestimation of modeled DO in spring even though the modeled temperature is underestimated, too. In this study, the biases between the observed and modeled results were corrected using a procedure adapted by Yara et al. (2011). Hence, Frontiers in Marine Science | www.frontiersin.org where f is the future modeled value of a parameter, and obs and p are the present observed and modeled values of the parameter, respectively.

Physical Parameters
Large seasonal fluctuation in the water temperature was found ( Figure 4A). The maximum and minimum temperatures during the monitoring period were 24.3 • C in August 2018 and 0.9 • C in February 2019, respectively. These were similar to observations at three sites at depths of 0 and 5 m in the same bay in the 1990s (Nakata et al., 2001). The salinity also fluctuated seasonally, but in a more complicated manner (Figure 4B). The maximum and minimum salinity during the monitoring period was 35.9 in July 2017 and 10.8 in November 2017, respectively. The salinity was relatively low from late winter to early spring as a result of the inflow of snow-melt water discharged directly to the bay from the surrounding land, as also mentioned in Nakata et al. (2001). The low salinity that appeared sporadically in November 2017 was first found by the continuous monitoring. This sudden low salinity was thought to be associated with an extreme storm on November 11, 2017.

Biogeochemical Parameters
The DO also showed strong seasonal fluctuation, and was generally high in boreal winter (January -March) and low in boreal summer (July -September; Figure 4C). The maximum and minimum DO during the monitoring period was 19.8 mg l −1 in April 2018 and 7.2 mg l −1 in August 2017, respectively. The seasonal fluctuation was characterized by the difference in oxygen saturation with temperature, and was high in cold water and low in warm water. The DO concentration is negatively correlated with temperature year-round, with a strong annual mean negative correlation of R = -0.86 (Table 1 and Figure 5A).
The DO also showed clear daily fluctuation, unlike the physical parameters. A maximum amplitude of 9.5 mg l −1 was recorded in May 2018 during the spring bloom ( Figure 4C; Yamaka, 2019). The maximum amplitude was comparable to or larger than reported in other coastal regions (e.g., 220 µmol kg −1 or around 6.9 mg l −1 ; Frieder et al., 2012). DO is generally higher in the daytime and lower at night (Figure 6A). The daily fluctuation in DO concentration is consistent with DO produced by photosynthesis by primary producers, such as phytoplankton, seaweeds, and seagrasses in the daytime and DO consumption by the respiration of all marine organisms at night. Therefore, the daily amplitude of DO concentration was relatively large in boreal spring (April -June) when photosynthesis was active, and the monthly-mean value was the largest in April (1.6 mg l −1 ; Figure 7A).
Accordingly, pH showed similar seasonal and daily fluctuations (Figures 4D,E, 6B,C). The maximum and minimum  (Nakata et al., 2001). Open dots in (D-F) denote in-situ values obtained by water samples. The dotted and solid lines in (C) denote the critical level for deoxygenation in Turner et al. (2005) and Vaquer-Sunyer and Duarte (2008; 2.0 mg l −1 ), and Marumo and Yokota (2012;4.29 mg l −1 ), respectively. Yellow and red domains in (f) denote the marginal and critical level of ocean acidification for calcifiers compiled in previous studies (between 1.1 and 1.5 and below 1.1 for arag , respectively; e.g., Waldbusser et al., 2015;Onitsuka et al., 2018). Bold values are less than -0.7 or larger than 0.7 (i.e., strong correlations). N/A: no data.
pH during the monitoring period was 8.72 (pH_25 • C = 8.35) in February 2019 and 7.86 (pH_25 • C = 7.71) in June 2018, respectively. Similar to DO, the seasonal and daily fluctuations were primarily characterized by physical and biogeochemical processes, respectively. The pH is dependent on temperature and is high in cold water and low in warm water, which contributes to the seasonal fluctuation in pH. The monitored pH_25 • C had a positive correlation with the temperature (an annual-mean correlation coefficient of R = 0.60; Table 1 and Figure 5B), although the correlation varied seasonally. Similar to the DO concentration, there was a daily fluctuation in pH, which was high in the daytime and low at night, due to photosynthesis by primary producers in the daytime and respiration by marine organisms at night (Figures 6B,C). The daily amplitude was relatively large in April, when photosynthesis is predominant, and the monthly-mean was the largest in April (0.11; Figure 7B), similar to that of the DO concentration. The pH also fluctuated sporadically. The maximum diurnal fluctuation in the study period was 0.64 (0.42 for pH_25 • C), which occurred in January 2019 (Figures 6C, 7C), and was comparable to fluctuations in estuarine (0.99), near-shore (0.49), and kelp forest (0.36-0.54) sites in other coastal regions (Hofmann et al., 2011;Frieder et al., 2012). These results indicate that coastal marine species in Oshoro Bay temporarily   experienced low pH resulting from the wide diurnal variation in pH, similar to other coastal regions.
The arag , another ocean acidification metric, also showed significant seasonal and daily fluctuations, but in a more complicated way (Figure 4F). The maximum and minimum arag was 4.7 in August 2017 and 1.3 in November 2017, respectively. The extremely low arag in November 2017 is thought to be associated with an extreme storm on November 11, 2017, which caused low temperature and salinity. The sporadic lowering of arag , which is close to the critical level of ocean acidification ( arag = 1.1; Onitsuka et al., 2018) was first found by multi-year monitoring. No arag values were available when either the salinity or pH data were not obtained. Therefore, it is likely that we failed to capture the annual minimum value of arag in this study. However, referring to observed seasonal fluctuations in arag (e.g., Ishii et al., 2011;Wakita et al., 2017Wakita et al., , 2021, the annual minimum should appear in winter when the water temperature is low and the value may have already reached close to the critical level of acidification. Similar to the DO concentration and pH, the daily fluctuation in arag was also prominent and it was high in the daytime and low at night, driven by daily fluctuation of carbon and oxygen dynamics through photosynthesis by primary producers in the daytime and respiration by marine organisms at night (Figure 6D). Although we did not cover the overall temporal fluctuation of arag because of missing data necessary for calculating the value (Figure 4F), the observed maximum daily fluctuation during the study period was 2.3 in April 2018, when the spring bloom was dominant. The monthlymean daily amplitude was as large as 0.6 in spring and summer ( Figure 7D).

Future Projections
The annual-mean projected rise in water temperature from the present to the end of this century was 2.2 • C for the RCP 2.6 scenario and 5.5 • C for the RCP 8.5 scenario ( Figure 8A). The rise in temperature was followed by significant lowering of the DO concentration, pH, and arag in the future (Figures 8B-D). However, the extent differed between the RCP scenarios: compared to their present states both physical and biogeochemical parameters are anticipated to change significantly for RCP 8.5, but much more moderately for RCP 2.6. The estimated extent of the projected results for RCP 2.6 is closer to the present values than those for RCP 8.5, which concurred with other future projections (e.g., Takao et al., 2015).
The projected annual-mean DO concentration decreased by 0.43 mg l −1 for RCP 2.6 and 1.02 mg l −1 for RCP 8.5, from the present (Figure 8B). Similarly, the projected annual-mean pH decreased by 0.06 for RCP 2.6 and 0.38 for RCP 8.5, from the present (Figure 8C). The projected pH reduction for RCP 2.6 and RCP 8.5 is comparable to that projected globally in previous studies (0.06-0.07 and 0.30-0.33, respectively; e.g.,  Hoshikawa, 2006;Shibano et al., 2014), both of which are main fisheries targets of the subarctic coast. Yellow and red domains in (D) denote the marginal and critical level of ocean acidification for calcifiers compiled in previous studies (between 1.1 and 1.5 and below 1.1 for arag , respectively; e.g., Waldbusser et al., 2015;Onitsuka et al., 2018). IPCC AR5, Jiang et al., 2019). The projected annual-mean arag decreased by 0.1 for RCP, 2.6 and 1.1 for RCP 8.5, from the present ( Figure 8D). The monitoring showed that the arag did not reach the critical level in the study period, but was sometimes marginal (Figure 4F). For the RCP 8.5 scenario, the monthly-mean arag is very often marginal in the 2090s. Considering that the future daily fluctuations of the biogeochemical parameters presumably magnify as the buffer capacity weakens, the daily minimum values will more often reach critical levels for calcifiers (<1.1 for arg ). Previous experimental results showed that the growth and malformation rates of larval Ezo abalone decreased and increased greatly, respectively, when the larvae were exposed in experimental seawater in which the arag was about 1.1, even for only 2-3 days (Onitsuka et al., 2018). This implies that more frequent future critical levels may severely damage calcifiers, especially during their early stages when they are relatively vulnerable to high CO 2 . This will also affect local industries, especially shellfisheries. Because around 60% of the shellfish and half of the calcifiers in Japan are caught in subarctic seas in which the arag is relatively low, ocean acidification is of great concern and the social demand for adaptive fisheries, such as changing the aquaculture style as mentioned in section "Suggestion for future adaptation, " in response to global change is high.

Criteria for Ocean Warming, Acidification, and Deoxygenation
Similarly, several studies have reported various critical levels for deoxygenation in which the DO concentration is <2.0 mg l −1 (Turner et al., 2005;Vaquer-Sunyer and Duarte, 2008) or 4.29 mg l −1 (Marumo and Yokota, 2012). Our monitoring results show that the minimum DO concentration during the study period (7.2 mg l −1 in August 2017) was still above the critical level ( Figure 4C). However, similar to pH and arag , if we assume that the daily fluctuation in DO concentration in the future is similar to that at present, i.e., the maximum daily fluctuation is as large as 9.5 mg l −1 , the daily minimum DO concentration might often reach critical levels of less than 2.0 mg l −1 in the future.
Along with ocean acidification and deoxygenation, concurrent ocean warming may also affect marine organisms, including calcifiers. For example, the main fisheries targets in the study site are vulnerable to warm temperatures. A representative critical temperature level was reported to be 23 • C for sea urchins and scallops (Hoshikawa, 2006;Shibano et al., 2014). Temperatures over 23 • C are not common at the study site now, but are projected to appear for 2 months for the RCP 2.6 scenario and 3 months for the RCP 8.5 scenario in summer ( Figure 8A). This means that calcifiers along the subarctic coast will suffer from ocean acidification in winter and ocean warming in summer in the future.

Suggestion for Future Adaptation
The large difference in the future projections between the RCP 2.6 and 8.5 scenarios clearly shows the need to reduce anthropogenic CO 2 and other greenhouse gas emissions as per the Paris Agreement, which has a target similar to the RCP 2.6 pathway, to alleviate fatal impacts of global change on marine organisms along subarctic coasts. It is also important to adapt to the upcoming impacts of global change. Adaptation strategies include modifying local industries in subarctic coastal regions, such as fisheries, and tourism such as eating around seafoods and recreational scuba diving. Because of overused fisheries resources, aquaculture plays a greater role in maintaining a steady supply of fisheries resources. Shellfish such as scallops, oysters, and abalone, and sea urchins are all calcifiers, and therefore, are concerned to be more or less affected by ocean acidification. However, the vulnerability to low pH and arag differs with various aspects such as the species abundance, life stage and habitats. For example, possible adaptation strategies to upcoming ocean acidification include changing the aquaculture style, such as raising larval calcifiers in low CO 2 conditions manipulated artificially, such as on-land aquaculture (e.g., Fujii, 2020), following the scientific insights that early larvae are relatively vulnerable to high CO 2 (e.g., Kurihara, 2008;Kimura et al., 2011;Waldbusser et al., 2015;Onitsuka et al., 2018;Foo et al., 2020). Combined impacts of global warming and ocean acidification are generally concerned to magnify the adverse effects on calcifiers, while global warming may partly and temporally alleviate the impacts of ocean acidification by rise in alkalinity from the land, as already evidenced in some other subarctic regions (e.g., Watanabe et al., 2009;Salisbury and Jönsson, 2018).
Regarding deoxygenation, low oxygen conditions are due to local human impacts, as well as global change and subsequent ocean warming. As local human impacts such as nutrification from the land are easier to alleviate than global change, further regulation of nutrient discharge to the coasts should delay the serious impacts of deoxygenation on marine organisms (e.g., Rabalais et al., 2010). However, the extent of regulation of nutrient discharge has not yet been quantified and further analytical studies are required.

CONCLUSION
This study was the first multi-year in-situ monitoring study of a subarctic coast and produced future projections of physical and biogeochemical parameters related to global change-driven ocean warming, acidification, and deoxygenation. The monitoring identified new characteristics of physical and biogeochemical parameters at the study site, such as sporadic changes in salinity and DO, which revealed that marine creatures sometimes experience marginal levels of ocean acidification. Our future projection results suggest that in the future high CO 2 world, calcifiers that are accustomed to current subarctic environments may be seriously damaged by ocean warming in summer and ocean acidification in mid-winter, and possibly by deoxygenation assuming that the daily fluctuation in DO concentration in the future is similar to that at present. Therefore, the mitigation of global change by following the Paris Agreement is of great importance for saving marine ecosystems and local industries, although local adaptation strategies are also required, such as adjusting aquaculture styles.

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