A Lacustrine Biomarker Record From Rebun Island Reveals a Warm Summer Climate in Northern Japan During the Early Middle Holocene Due to a Stronger North Pacific High

The summer climate of northern Japan since the last glacial period has likely been determined by atmospheric and oceanic dynamics, such as changes in the North Pacific High, the position of the westerlies, the Kuroshio Current, the Tsushima Warm Current (TWC), and the East Asian summer monsoon. However, it is unclear which factor has been most important. In this study, we analyzed leaf wax δ13C and δD and glycerol dialkyl glycerol tetraethers (GDGTs) in sediments from Lake Kushu, Rebun Island, northern Japan, and discuss changes in climate over the past 17,000 years. The GDGT-based temperature, the averaged chain length, δ13C and δD of long-chain n-fatty acids indicated that the climate was cold during the Oldest Dryas period ∼16 ka and warm in the early Middle Holocene from ∼9 to 6 ka. This climate change is consistent with the sea surface temperature in the Kuroshio–Oyashio transition, but inconsistent with changes in the TWC in the Sea of Japan. The results imply that the summer climate of northern Japan was controlled mainly by changes in the development of the North Pacific High via changes in the position of the westerly jet and East Asian summer monsoon rainfall, whereas the influence of the TWC was limited over a millennial timescale.


INTRODUCTION
The summer climate of northern Japan is influenced by the summer position of the westerly winds and the strength of the Tsushima Warm Current (TWC). The position of the westerly winds is linked to the oceanic subarctic boundary between the subtropical Kuroshio and subarctic Oyashio currents in the Pacific (Figure 1) (e.g., Yamamoto et al., 2004;Yamamoto et al., 2005;Isono et al., 2009;Yamamoto, 2009). Yamamoto et al. (2005) and Isono et al. (2009) identified the southern summer position of the Kuroshio-Oyashio transition (KOT) in the Oldest Dryas and Younger Dryas periods due to the weaker North Pacific High, a northward shift of the KOT at the end of the Younger Dryas period, the northernmost position at 8 ka and the gradual southward shift during the Middle and Late Holocene. The TWC flows northwards as a branch of the Kuroshio Current along the eastern margin of the Sea of Japan (Figure 1). The TWC promotes moisture uptake by the predominant winter monsoon winds and carries heat to northern Japan. Diatom records from four offshore sites along the Sea of Japan coast showed that TWC species appeared after 7 ka and fluctuated in abundance on a millennial timescale (Koizumi et al., 2006).
Pollen records from Hokkaido (Igarashi, 2013 and references therein) and Lake Kushu (Müller et al., 2016;Leipe et al., 2018) on Rebun Island showed that the vegetation of northern Japan changed from a boreal forest in the Late Glacial to a cooltemperate forest in the Early Holocene due to deglacial warming. Temperate deciduous oak trees were more widely spread during the Middle Holocene, while forests containing more Abies and Pinaceae trees spread during the Late Holocene due to gradual cooling.
Other proxies such as Sphagnum and vascular plant cellulose δ 18 O in peat are more sensitive to summer hydrological conditions. They showed remarkable changes in the high moors in Rishiri (Sites MHWL-1 and -3) and Hokkaido (Site BKB-2) during the Late Holocene, which corresponded to changes in the summer sea surface temperature (SST) in the KOT and the TWC strength (Sakurai et al., 2021;Yamamoto and Seki, unpublished data).
The summer westerly wind axis, moisture content, temperature, precipitation, and precipitated water δD and d-excess seasonally covaried at the study site (Figure 2), indicating a close relationship between the position of the summer westerly wind axis, the East Asian summer monsoon, and precipitated water δD and d-excess ( Figure 2). In summer, the westerly wind axis shifts north, a warm and moist air mass occupies the study site (45°N), the δD of precipitated water is heavier, and the d-excess is lower. In winter, the westerly wind axis shifts south, a cold and dry air mass occupies the study site, the δD of precipitated water is lighter, and the d-excess is higher. The higher d-excess in winter is attributed to kinetic fractionation during the evaporation process in the Sea of Japan under a large humidity deficit and a strong temperature contrast in the air-sea surface interface (Sugimoto et al., 1988;Li et al., 2017).
In this study, we analyzed leaf wax δ 13 C and δD and glycerol dialkyl glycerol tetraethers (GDGTs) in sediments from Lake Kushu, Rebun Island, northern Japan, and discuss changes in the summer climate and the atmospheric and ocean dynamics that affected the climate of northern Japan over the past 17 ka.

Sediment and Chronology
Lake Kushu is located in the northern part of Rebun Island (45°25′58″N, 141°02′05″E) ( Figure 1), 230-400 m from the modern sea coast. The bean-shaped lake has a maximum length of ∼1,100 m. The maximum water depth is ∼6 m in the eastern part of the lake, with average depths of 3-5 m. In February 2012, when a thick ice layer covered the lake, Dokon (Sapporo, Japan) performed coring in the central part of the lake. Composite core RK12 was recovered from two boreholes (RK12-1 and RK12-2) located within a few meters of each other.
Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 704332 laminated with sections of relatively low to high organic matter concentrations. The upper 850 cm of the sediment column consists of homogeneous organic-rich clay.
An age-depth model of the RK12 core was established based on 57 accelerated mass spectrometry (AMS) radiocarbon dates (Müller et al., 2016), showing continued sedimentation beginning ∼17 ka  Table S1). All available AMS dates were converted into calendar ages using Oxcal v4.4 (Bronk Ramsey, 2021) and the latest IntCal20 calibration curve (Reimer et al., 2020). The 95.4% probability ranges obtained in this study (Supplementary Table S1) show very minor differences from the data published in Müller et al. (2016), thus indicating general robustness of the published age model. For the current study, an updated Bayesian age model for the Kushu sequence has been constructed using Oxcal v4.4 (Bronk Ramsey, 2021) based on the 57 AMS 14 C dates (Müller et al., 2016) and three tephra ages using Oxcal v4.4 (Bronk Ramsey, 2021). The model utilizes a P_Sequence deposition model (Bronk Ramsey, 2008), with a variable k parameter (Bronk Ramsey and Lee, 2013) and a "General" Outlier_Model (Bronk Ramsey, 2009), applying the latest IntCal20 calibration curve (Reimer et al., 2020). The 14 C dates regarded as outliers (see Müller et al., 2016 and Supplementary Table S1) have been removed from the model. Three precise tephra ages have been incorporated into the model for providing additional chronological constraints, given the identification of these tephra layers in the RK12 sediment core of Lake Kushu (Chen et al., 2016;Chen et al., 2019). These include B-Tm tephra (composite core depth-150.5 cm; 1,004 cal yr BP; Oppenheimer et al., 2017), Ko-g tephra (core depth-1,169 cm; 95.4% probability range-6,651-6,446 cal yr BP; Chen et al., 2021) and Ma-f∼j tephra (core depth-1,277 cm; 95.4% probability FIGURE 3 | Distribution of 57 calibrated AMS radiocarbon and three tephra dates used to construct Bayesian P_sequence depositional age-depth model for the RK12 sediment core from Lake Kushu (Supplementary Table S1).
In the upper half of the core, representing the past ∼6 ka, the linear sedimentation rate is very high (i.e., about 1 cm every 6 years). In the bottom half of the core, the age-depth model shows greater variation in the sedimentation rates and substantially slower sedimentation (i.e., about 1 cm every 20 years) before ∼9.5 ka. Based on a multi-proxy approach (diatoms, aquatic pollen, algae, geochemistry), Schmidt et al. (2019) reconstructed three phases of lake basin development, including a marshy phase between 16.6 and 9.4 cal ka BP, a lagoon phase between 9.4 and 5.9 cal ka BP, and a freshwater lake phase after 5.9 cal ka BP, as the sand bar separating the Kushu lagoon from the sea formed.
The concentration of n-fatty acids in the acid fraction was analyzed using gas chromatography (GC) with an Agilent 6890 series gas chromatograph with on-column injection and electronic pressure control systems, and a flame ionization detector (Agilent, Santa Clara, CA, United States). Samples were dissolved in toluene. Helium was the carrier gas and the flow velocity was maintained at 30 cm/s. An Agilent J&W CP-Sil 5 CB column was used (length 50 m; i.d. 0.25 mm; thickness 0.25 μm). The oven temperature was programmed from 100 to 130°C at 20°C/min and from 130 to 310°C at 4°C/min, and then maintained at 310°C for 30 min. The average chain length (ACL) of n-fatty acids in this study is defined as: The δ 13 C (‰ Vienna Pee Dee Belemnite [VPDB]) of the methylated n-fatty acids in the acid fraction was analyzed using a Thermo Fisher Scientific (Waltham, MA, United States) GC IsoLink II with a capillary column coated with a TG-5MS column (30 m length; i.d. 0.25 mm; 0.25 μm film thickness) combined with a Delta V mass spectrometer through a combustion furnace at 1,000°C. The sample was dissolved in toluene. The oven temperature was programmed from 100 to 230°C at 30°C/min and from 230 to 310°C at 4°C/min and then maintained at 310°C for 15 min. As an internal isotopic standard, n-C 36 H 74 was used to check the measurement conditions. Data were converted into values relative to VPDB using standard delta notation by comparison with CO 2 standard gas. The δ 13 C values of fatty acids were obtained from the measured values of fatty acid-methyl esters by correcting for methyl carbon (−34.1‰). The reproducibility of the measurement based on repeated analyses was better than ±0.1‰.
The neutral fraction was separated into four fractions using SiO 2 column chromatography: F1, 3 ml hexane; F2, 3 ml hexane:toluene (3:1); F3, 4 ml toluene; and F4, 3 ml toluene:CH 3 OH (3:1). An aliquot of F4 was dissolved in hexane-2-propanol (99:1), spiked with an internal standard (500 ng of C 46 glycerol trialkyl glycerol tetraether [GTGT]; Patwardhan and Thompson, 1999), and filtered through a short bed of Na 2 SO 4 . GDGTs were analyzed using highperformance liquid chromatography-mass spectrometry (HPLC-MS) with an Agilent 1260 HPLC system coupled to an Agilent 6130 Series quadrupole mass spectrometer. Separation was achieved on two ultra-HPLC silica columns (ACQUITY BEH HILIC columns, 2.1 × 150 mm, 1.7 μm; Waters, Milford, MA, United States) in series, fitted with a 2.1 × 5 m pre-column of the same material (Waters) and maintained at 30°C. GDGTs were eluted isocratically for 25 min with 18% B, followed by a linear gradient to 35% B in 25 min, then a linear gradient to 100% B in 30 min, where A was hexane and B was hexane:isopropanol (9:1, v/v). The flow rate was 0.2 ml/min. The total run time was 90 min with a 20 min re-equilibration. Ionization was achieved using atmospheric pressure, positive ion chemical ionization-MS. The spectrometer was run in selected ion monitoring mode (m/z 743.8, 1,018, 1,020, 1,022, 1,032, 1,034, 1,036, 1,046, 1,048, 1,050, 1,292.3, 1,296.3, 1,298.3, 1,300.3, and 1,302.3). Compounds were identified by comparing mass spectra and retention times with those in the literature (Hopmans et al., 2000;De Jonge et al., 2014). Quantification was achieved by integrating the peak area in the (M+H) + chromatogram and comparing these with the peak area of an internal standard (C 46 GTGT) in the (M+H) + chromatogram according to the method of Huguet et al. (2006). The correction value of ionization efficiency between GDGTs and the internal standard was obtained in our laboratory by comparing the peak areas of isolated crenarchaeol and branched GDGTs I and II (Schouten et al., 2013) and C 46 GTGT in known amounts. In the routine analysis, a working standard that was a mixture of C 46 GTGT and the GDGTs extracted and purified from East China Sea sediment was inserted every 20 samples to monitor ionization efficiency changes.
To estimate the terrestrial/in-situ production ratio of branched GDGTs, the weighted average number of cyclopentane moieties of tetramethylated branched GDGTs (

GDGT-Based Temperature
The MAT reconstruction based on the four different calibrations yielded similar variation within the calibration error of 4.6-5.0°C (Figure 4). The average value of the reconstructed temperature varied between -0.7 and 10.6°C and was low at ∼17 ka and high between 9 and 6 ka ( Figure 4). The core-top temperature (3.3-5.8°C) was consistent with the current mean annual air temperature ∼83 km southeast of the study site (6.4°C at Teshio; Japan Meteorological Agency, https://weather.time-j. net/Stations/JP/teshio) within the calibration error. Because branched GDGTs are produced in both terrestrial and FIGURE 5 | Abundance of pollen in sediments from the RK12 core during the last 17 ka (Müller et al., 2016).
Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 704332 7 lacustrine environments (e.g., Schouten et al., 2012), the reconstructed temperature may reflect both terrestrial and lake water temperatures. The #Rings tetra index varied between 0.1 and 0.5, but was mainly 0.3-0.4 after 16 ka (Figure 4), corresponding to the range of dominant terrestrial branched GDGTs (Sinninghe Damsté, 2016). This indicated that the branched GDGTs were mainly of terrestrial soil origin, with a small lacustrine contribution after 16 ka. We assume that the GDGT-derived MAT reflects the terrestrial air temperature rather than the lake water temperature. The record shows a colder climate in the Late Glacial and a warmer climate in the early Middle Holocene (Figure 4).
The pollen assemblage in sediments from the RK12 core was described and the paleovegetation discussed in Müller et al. (2016) and Leipe et al. (2018) is as follows. From 16 to 11 ka, the vegetation shifted from herbaceous (predominantly Cyperaceae) to trees and shrubs (Pinaceae, Alnus, and Betula) ( Figure 5). The spread of woody vegetation was interrupted during the Younger Dryas period, although the reversal was not prominent. In the Middle Holocene between 8 and 4 ka, there was a major spread of deciduous Quercus trees, with the highest percentages of Quercus pollen (17-27%) between 6 and 5 ka ( Figure 5) (Müller et al., 2016;Leipe et al., 2018). This is in line with a major spread of cool mixed and cool conifer forests in the Hokkaido Region (Igarashi, 2013). The Late Holocene section of the RK12 record showed the highest percentages of Abies and Pinaceae pollen by ∼2 ka (Figure 5), suggesting an increase in coniferous tree cover in the regional forest vegetation. This pollen record in Lake Kushu is consistent with the result of our MAT reconstruction. The culmination of Quercus pollen was delayed behind the MAT maximum by ∼1 ka, presumably reflecting the time required for vegetation transition.

Distribution of n-Fatty Acids
All of the samples showed a bimodal homologous distribution of n-fatty acids with maxima at C 16 and C 26 or C 28 . The n-fatty acids around C 16 are ubiquitous in eukaryotes and bacteria, while the longchain n-fatty acids (>C 26 ) are derived from vascular plants (Eglinton and Hamilton, 1967). The long-chain n-fatty acids (>C 26 ) were more abundant than the shorter homologs. The even carbon number preference ranged from 5 to 18, typical fresh leaf wax.
The ACL varied between 27.5 and 28.2 and was higher in the period between 10 and 6 ka than in other periods (Figure 4). In modern plants, the ACL of plant leaf wax (long-chain n-alkanes) is higher in warmer climates (Tipple and Pagani, 2013). It is also sensitive to relative humidity, but the response is opposite in different species (Hoffmann et al., 2013). The ACL is higher in grasses than trees in the temperate climate zone (Bush McInerney, 2013). The abundance variation in Poaceae and Cyperaceae pollen was not consistent with the ACL (Figure 5; Müller et al., 2016), suggesting that the higher ACL between 10 and 6 ka does not reflect grass abundance and aridity but rather a warm climate, as indicated by the GDGT indices.
Relatively high ACL was also found in the Oldest Dryas period (Figure 4). The pollen record indicated that Cyperaceae pollen was abundant in this period ( Figure 5; Müller et al., 2016), suggesting that the high ACL reflected an abundance of sedge.

Leaf Wax δ 13 C
The δ 13 C of long-chain (C 26 , C 28 , and C 30 ) n-fatty acids varied between -35.6 and -26.8‰ (Figure 4). The values were higher in the period between 9 and 5 ka (Figure 4). The δ 13 C of long-chain n-fatty acids of C 3 and C 4 plants sampled in the early 2000s (CO 2 concentration ∼370 ppm) averaged -37.1 ± 2.0‰ (n 13) and -19.5 ± 1.8‰ (n 9), respectively (Chikaraishi et al., 2004). However, these end-member values have varied in the past. Contemporary δ 13 C values for both C 3 and C 4 plants have been affected by the Suess effect (decreasing by ∼2‰ over the past 250 years; Keeling et al., 2017). The δ 13 C values of C 3 plants have been reduced further by increasing isotopic fractionation, which is governed by the atmospheric CO 2 concentration (Schubert and Jahren, 2015). Considering these effects, the calculated fatty acid δ 13 C of C 3 plants at a CO 2 of 280 ppm is -33.7‰. In comparison, the δ 13 C value for C 4 plants is independent of the atmospheric CO 2 concentration. The fatty acid δ 13 C of C 4 plants before 1750 CE was -17.5‰. The δ 13 C of long-chain n-fatty acids in the study samples varied around the end-member values of C 3 plants. If the variation in δ 13 C reflected the mixing ratio of C 3 and C 4 plants, the higher δ 13 C between 10 and 5 ka would indicate the contribution of C 4 plants. However, the pollen record from the study core did not show evidence of a greater contribution of C 4 plants in this period, suggesting that another factor was responsible for the reconstructed shift in the δ 13 C values.
The δ 13 C of C 3 plants is affected by photosynthetic activity (Farquhar and Richards, 1984;Evans et al., 1986), where greater activity induces higher δ 13 C. In the RK12 core, higher δ 13 C appeared in a warm period (Figure 4). Therefore, we postulate that a warm climate enhanced C 3 plant photosynthesis, which increased δ 13 C.

Leaf Wax δD
The δD of long-chain (C 26 , C 28 , and C 30 ) n-fatty acids varied between -220.6 and -151.5‰ (Figure 4). The δD was lower at ∼17, 11, and 5 ka and higher at 9-6 ka ( Figure 4). The variation in δD was similar to that in the GDGT-based MATs (Figure 4). Dansgaard (1964) and Jouzel et al. (1994) showed that the δ 18 O of global precipitation is related to the MAT. Jouzel et al. (1994) also found that areas with a MAT below 15°C had a linear relationship between MAT and the mean annual δ 18 O of precipitation (δ 18 O 0.64 ). This is equivalent to the formula δD 5.12 . The MAT at the study site is 6.4°C (1981-2010; Japan Meteorological Agency, https://weather.time-j.net/ Stations/JP/teshio) and the mean annual δD at Teshio 80 km southeast of the sampling site is −69.6‰ (2010-2014Li et al., 2017). These values are consistent with the relationship described above.
The maximum range of long-chain fatty acid δD variation was 69.1‰ (Figure 4). If this variation was caused only by a process following the above MAT-δD relationship, a 69.1‰ variation in δD would correspond to a 13.5°C variation in MAT. Our MAT estimates based on the GDGT compositions indicated that the temperature variation over the last 17 ka was 11.3°C. Summer SST records from the northern Japan margin over the past 17 ka show a variation of ∼8°C (Site PC-6 in Minoshima et al., 2007;Site GH02-1030in Inagaki et al., 2009. This implies that the temperature effect explains a large part of the variation observed in C 26 and C 28 fatty acid δD. Two other potential factors control long-chain fatty acid δD via the mean annual δD of precipitation. The first factor is the seasonal variation in the δD of precipitation. Water precipitated in autumn has ∼30‰ higher δD values than that precipitated in winter and spring (Figure 2) (Li et al., 2017). Thus, a large amount of autumn precipitation or a lower amount of winter and spring precipitation would shift the mean annual δD values in a positive direction.
The position of the westerly wind axis is the second factor controlling long-chain fatty acid δD via the mean annual δD of precipitation (Li et al., 2017). Precipitation occurs in Hokkaido when a low-pressure system migrates from west to east along with the westerly jet. The δD values of precipitated water were 24-32‰ higher when the migration path of the low-pressure system was located to the north of the study location relative to times when the migration path was located to the south (Li et al., 2017). When the westerly jet is located north of the study site, southerly winds bring moisture from the south, inducing more precipitation at the study site ( Figure 2). Thus, a FIGURE 6 | Changes in the glycerol dialkyl glycerol tetraether-based mean annual temperatures (MATs) obtained using different calibrations (averages and 1σ range), U K 37 '-based sea surface temperature (SST) at sites MD01-2421 (Yamamoto et al., 2004;Yamamoto et al., 2005;Isono et al., 2009;Yamamoto, 2009) and PC-6 (Minoshima et al., 2007) in the Kuroshio-Oyashio transition zone, temperature index of diatoms (Td') and relative abundance of Fragilariopsis doliola, and Tsushima Warm Current species in cores KH-84-3-9, KH-84-3-33, KH-86-2-9, and D-GC-6 (Koizumi et al., 2006) over the last 17 ka. northward shift in the position of the westerly jet increases the δD of precipitated water, and the associated higher fraction of summer-to-annual precipitation further increases the δD. These results suggest that the position of the westerly jet controls the amount of precipitation, as well as the δD of precipitated water over a long timescale. Higher long-chain fatty acid δD reflects a northerly position of the summer westerly jet and the resultant increase in summer and autumn precipitation.
One may speculate that the strengthening of the TWC increased winter precipitation and decreased the fatty acid δD from the Middle to Late Holocene because, in winter, the East Asian winter monsoon influences all of Hokkaido and brings heavy snowfall to mountainous areas. Plants might partly use the water supplied by winter precipitation. Even if winter precipitation was higher, the δD should be lower. However, the precipitation from January to April is only one third that from June to October at present (Japan Meteorological Agency, https://weather.time-j.net/Stations/JP/ teshio). If winter precipitation was lower in the Middle Holocene than in the Late Holocene, the effect on the δD of annual precipitation should be small. Therefore, this case can likely be neglected on a millennial timescale.
Evapotranspiration induces the enrichment of deuterium in the plant body water. The evapotranspiration rate depends on the relative humidity. When the relative humidity is low, evapotranspiration is expected to be more active, resulting in higher leaf wax δD. However, the low abundance of Poaceae and Cyperaceae pollen does not indicate a dry environment during 9-6 ka ( Figure 5; Müller et al., 2016) and was not consistent with the high δD during 9-6 ka, suggesting that the higher high δD does not reflect a dry environment.
The above discussion summarizes that the leaf wax δD at the study site reflects the global distribution of precipitated water δD and temperature. The latitudinal gradient of precipitated water δD is enhanced north and south of the westerly jet (Li et al., 2017). We thus conclude that the higher leaf wax δD during 9-6 ka likely reflected a warmer climate and the northern position of the westerlies.

Factors Controlling Summer Climate in Northern Japan Over the Past 17 ka
The Hokkaido region is located at the mean position of the modern summer westerly jet ( Figure 1) and at the northern margin of the East Asian summer monsoon, which brings a warm, moist air mass from the Pacific Ocean to this area. When the westerly jet is located in the north, the southerly winds (i.e., the East Asian summer monsoon) reach Hokkaido, resulting in a warm, wet climate (Nitta, 1987;Kawamura et al., 1998).
The GDGT-based MATs, the ACL, and the δ 13 C and δD of long-chain n-fatty acids indicated that climate was cold in the Oldest Dryas period and warm during 10-6 ka ( Figure 4). This climate change is consistent with the SST changes in the KOT ( Figure 6) (Site PC-6 in Minoshima et al. (2007); Site MD01-2421 in Yamamoto et al. (2004), Yamamoto et al. (2005), Isono et al. (2009), Yamamoto, (2009). The summer position of the westerly wind axis in the study region is related to the strength of the North Pacific High (Nitta, 1987;Kawamura et al., 1998), with a stronger (weaker) North Pacific High shifting the westerly wind axis northward (southward) (Figure 7). The development of the North Pacific High drives the oceanic subtropical gyre circulation in the North Pacific. The Pacific SST in the KOT is a good indicator of the strength of the Kuroshio and its extension, which form part of the subtropical gyre circulation. Thus, the correspondence between the climate of Rebun Island and the SST in the KOT suggests that the summer position of the westerly wind axis is linked to the oceanic subtropical gyre circulation in the North Pacific via regional atmospheric dynamics in the northwestern Pacific region. The postglacial climate of Rebun Island reflects changes in the atmospheric conditions, i.e., the northwardshifted position of the westerly jet, stronger influence of the FIGURE 7 | Schematic maps show the postulated positions of the summer westerly jet, the Kuroshio Extension (Yamamoto, 2009) and the Tsushima Warm Current (TWC; Koizumi et al., 2006) in different climate regimes.
Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 704332 East Asian summer monsoon, and development of the North Pacific High, and vice versa (Sakurai et al., 2021). However, the climate change during the Holocene is not consistent with changes in the TWC in the Sea of Japan indicated by the temperature index of diatoms [Td' warm/ (warm + cold) water diatom species] or the abundance of Fragilariopsis doliola, a marker diatom species of the TWC, in core D-GC-6 from the southern Sea of Japan ( Figure 6) (Koizumi et al., 2006). On Rishiri Island, 40 km south of Rebun Island, cellulose δ 18 O in peat cores from site MHWL showed that the influence of the TWC on the δ 18 O of precipitated water appeared after 2.8-1.3 ka (Yamamoto and Seki, unpublished data). However, the relationship between the fatty acid δD and the strength of the TWC on a millennial timescale during the Holocene is not clear.
The above considerations imply that the climate of northern Japan was controlled mainly by changes in the development of the North Pacific High via changes in the position of the westerly jet and the East Asian summer monsoon rainfall, while the influence of the TWC was limited on a millennial timescale.

CONCLUSION
The GDGT-based MAT, ACL, δ 13 C and δD of long-chain n-fatty acids indicated that the climate in northern Japan was cold in the Oldest Dryas period and warm during 9-6 ka. This climate change is consistent with the SST in the KOT, but inconsistent with changes in the TWC in the Sea of Japan. The results imply that the climate of northern Japan was controlled mainly by changes in the development of the North Pacific High via changes in the position of the westerly jet and the East Asian summer monsoon rainfall, whereas the influence of the TWC was limited on a millennial timescale.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
MY conceived the research idea, contributed biomarker data and prepared the manuscript. FW and KS contributed biomarker data. KY and TH conducted the field work. KG, HY and PT managed the drilling project. PT and HYC updated the age model. All authors discussed the data and improved the manuscript.

FUNDING
This study was supported by the Japan Society for the Promotion of Science (25610146, JPMXS05R2900001, and 19H05595 to MY and 26101002, 21101002 to HY).