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 δ¹³C 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, δ¹³C 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).

FIGURE 1

FIGURE 1. Map showing study site RK12-2 and the reference site locations, i.e., MHWL-1 and 3 (Yamamoto and Seki, in preparation), BKB-2 (Sakurai et al., 2021), PC-6 (Minoshima et al., 2007), MD01-2421 (Yamamoto et al., 2004), KH-84-3-9, KH-84-3-33, KH-86-2-9, and D-GC-6 (Koizumi et al., 2006), and D-GC-9 (Koizumi et al., 2006), as well as the mean positions of the summer (June to August) westerly wind axis with the range of variation for 30 years from 1988 to 2017 (yellow shading) of summer westerlies, the Kuroshio, Kuroshio Extension, Oyashio, Tsushima Warm Current (TWC), and Tsugaru Warm Current (T). T indicates the meteorological station at Teshio.

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 cool-temperate 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 δ¹⁸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).

FIGURE 2

FIGURE 2. (A) Zonal wind velocity at 200 hPa (positive westward) along 145°E, (B) Temperature and moisture content at 850 hPa along 145°E (Sakurai et al., 2021), (C) Precipitated water δD, d-excess, temperature and precipitation at the site of Teshio Town (2010–2013; Li et al., 2017) 80 km southeast of the study site, and the position of the westerly wind axis. Black solid dots indicate the positions of the westerly wind axis.

In this study, we analyzed leaf wax δ¹³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.

Samples and Methods

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.

The composite core revealed a continuous, partly laminated, organic-rich 20.3 m-long sediment column. The lithology is described in Müller et al. (2016). In brief, the basal unit contains sandy clay with pebbles, suggesting a stronger river influence. Peat (1,925–1,935 and 1,905–1,915 cm) and organic-rich clay (1,915–1,925 and 1,895–1,905 cm) layers appear in the lower part of the core, suggesting shallow water or marshy environments. The interval within 1,895–1,390 cm is characterized by homogeneous, relatively organic-poor clay, which is interrupted only by two sand seams (1,790–1,815 and 1,765–1,780 cm). Within 1,390–850 cm, the clay is mostly finely 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 (Supplementary 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 ¹⁴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 ¹⁴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 range—7,550–7,128 cal yr BP; Chen et al., 2021). The 95.4% Highest Probability Density (HPD) range for the deposition model is illustrated in Figure 3.

FIGURE 3

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.

Analytical Methods

Lipids were extracted (×3) from ∼1 g of dried sediment using a Dionex ASE 200 Accelerated Solvent Extractor (Dionex, Sunnyvale, CA, United States) with dichloromethane:methanol (6:4). The extract was separated into neutral and acid fractions via aminopropyl silica gel column chromatography (Gao et al., 2015). The acid fraction was methylated with methanol:acetyl chloride (95:5) and purified with SiO₂ column chromatography.

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:

$\text{ACL}=\frac{\left({26\text{C}}_{26}+{28\text{C}}_{28}+{30\text{C}}_{30}+{32\text{C}}_{32}\right)}{\left({\text{C}}_{26}.+{\text{C}}_{28}+{\text{C}}_{30}+{\text{C}}_{32}\right)}.$

The δ¹³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₃₆H₇₄ was used to check the measurement conditions. Data were converted into values relative to VPDB using standard delta notation by comparison with CO₂ standard gas. The δ¹³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₂ column chromatography: F1, 3 ml hexane; F2, 3 ml hexane:toluene (3:1); F3, 4 ml toluene; and F4, 3 ml toluene:CH₃OH (3:1). An aliquot of F4 was dissolved in hexane-2-propanol (99:1), spiked with an internal standard (500 ng of C₄₆ glycerol trialkyl glycerol tetraether [GTGT]; Patwardhan and Thompson, 1999), and filtered through a short bed of Na₂SO₄. GDGTs were analyzed using high-performance 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₄₆ 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₄₆ GTGT in known amounts. In the routine analysis, a working standard that was a mixture of C₄₆ 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 (#Rings_tetra) was calculated (Sinninghe Damsté, 2016):

$\#{\text{Rings}}_{\text{tetra}}=\frac{\left(\left[\text{Ib}\right]+2\left[\text{Ic}\right]\right)}{\left(\left[\text{Ia}\right]+\left[\text{Ib}\right]+\left[\text{Ic}\right]\right)}.$

To estimate the mean annual temperature (MAT), the following parameters were calculated with the following equations based on a global soil dataset obtained using different calibration strategies (De Jonge et al., 2014):

$\begin{array}{l}\text{MB}{T}_{5\text{ME}}^{\prime }=\frac{\left(\left[\text{Ia}\right]+\left[\text{Ib}\right]+\left[\text{Ic}\right]\right)}{\left(\left[\text{Ia}\right]+\left[\text{Ib}\right]+\left[\text{Ic}\right]+\left[\text{IIa}\right]+\left[\text{IIb}\right]+\left[\text{IIc}\right]+\left[\text{IIIa}\right]\right)}\text{,}\\ \text{MAT}=-8.57+31.45\text{MB}{T}_{5\text{ME}}^{\prime }\\ (\text{RMSE}=4.8\mathrm{°}\text{C})\end{array}(1)$

$\begin{array}{l}\text{Index}1=\log \frac{\left(\left[\text{Ia}\right]+\left[\text{Ib}\right]+\left[\text{Ic}\right]+\left[\text{II}{a}^{\prime }\right]+\left[\text{III}{a}^{\prime }\right]\right)}{\left(\left[\text{Ic}\right]+\left[\text{IIa}\right]+\left[\text{IIc}\right]+\left[\text{IIIa}\right]+\left[\text{III}{a}^{\prime }\right]\right)}\text{,}\\ \text{MAT }=5\text{.}05+14\text{.}86\text{Index }1\\ (\text{RMSE }=4\text{.}7\text{°C})\end{array}(2)$

$\begin{array}{l}\text{MATmr }=7\text{.}71+17\text{.}1\left[\text{Ia}\right]+25\text{.}9\left[\text{Ib}\right]+34\text{.}4\left[\text{Ic}\right]-28\text{.}6\left[\text{IIa}\right]\text{,}\\ \text{RMSE }=4\text{.}6\text{°C}\text{.}\end{array}(3)$

$\begin{array}{l}\text{MATmrs }=5\text{.}58+17\text{.}91\left[\text{Ia}\right]-18\text{.}77\left[\text{IIa}\right]\\ (\text{RMSE }=5\text{.}0\text{°C})\end{array}(4)$

where [Ia] to [IIIc′] are the GDGTs defined in De Jonge et al. (2014).

Results and Discussion

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 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).

FIGURE 4

FIGURE 4. Changes in the glycerol dialkyl glycerol tetraether (GDGT)-based mean annual temperatures (MATs) obtained using different calibrations (Equations 1–4), weighted average number of cyclopentane moieties of tetramethylated branched GDGTs (#Rings_tetra), average chain length (ACL), and δ¹³C and δD of long-chain n-fatty acids in sediments from core RK12-2, Kushu Lake, Rebun Island, over the last 17 ka. The vertical bar of the sample value indicates 1σ. OD, Oldest Dryas period; BA, Bølling–Allerød period; and YD, Younger Dryas period.

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.

FIGURE 5

FIGURE 5. Abundance of pollen in sediments from the RK12 core during the last 17 ka (Müller et al., 2016).

Distribution of n-Fatty Acids

All of the samples showed a bimodal homologous distribution of n-fatty acids with maxima at C₁₆ and C₂₆ or C₂₈. The n-fatty acids around C₁₆ are ubiquitous in eukaryotes and bacteria, while the long-chain n-fatty acids (>C₂₆) are derived from vascular plants (Eglinton and Hamilton, 1967). The long-chain n-fatty acids (>C₂₆) 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 δ¹³C

The δ¹³C of long-chain (C₂₆, C₂₈, and C₃₀) 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 δ¹³C of long-chain n-fatty acids of C₃ and C₄ plants sampled in the early 2000s (CO₂ 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 δ¹³C values for both C₃ and C₄ plants have been affected by the Suess effect (decreasing by ∼2‰ over the past 250 years; Keeling et al., 2017). The δ¹³C values of C₃ plants have been reduced further by increasing isotopic fractionation, which is governed by the atmospheric CO₂ concentration (Schubert and Jahren, 2015). Considering these effects, the calculated fatty acid δ¹³C of C₃ plants at a CO₂ of 280 ppm is –33.7‰. In comparison, the δ¹³C value for C₄ plants is independent of the atmospheric CO₂ concentration. The fatty acid δ¹³C of C₄ plants before 1750 CE was –17.5‰. The δ¹³C of long-chain n-fatty acids in the study samples varied around the end-member values of C₃ plants. If the variation in δ¹³C reflected the mixing ratio of C₃ and C₄ plants, the higher δ¹³C between 10 and 5 ka would indicate the contribution of C₄ plants. However, the pollen record from the study core did not show evidence of a greater contribution of C₄ plants in this period, suggesting that another factor was responsible for the reconstructed shift in the δ¹³C values.

The δ¹³C of C₃ plants is affected by photosynthetic activity (Farquhar and Richards, 1984; Evans et al., 1986), where greater activity induces higher δ¹³C. In the RK12 core, higher δ¹³C appeared in a warm period (Figure 4). Therefore, we postulate that a warm climate enhanced C₃ plant photosynthesis, which increased δ¹³C.

Leaf Wax δD

The δD of long-chain (C₂₆, C₂₈, and C₃₀) 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 δ¹⁸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 δ¹⁸O of precipitation (δ¹⁸O = 0.64 MAT—12.8). This is equivalent to the formula δD = 5.12 MAT—92.4. 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–2014; Li 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-1030 in Inagaki et al., 2009). This implies that the temperature effect explains a large part of the variation observed in C₂₆ and C₂₈ 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 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 δ¹³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 northward-shifted position of the westerly jet, stronger influence of the East Asian summer monsoon, and development of the North Pacific High, and vice versa (Sakurai et al., 2021).

FIGURE 6

FIGURE 6. Changes in the glycerol dialkyl glycerol tetraether-based mean annual temperatures (MATs) obtained using different calibrations (averages and 1σ range), U^K₃₇'-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.

FIGURE 7

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.

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 δ¹⁸O in peat cores from site MHWL showed that the influence of the TWC on the δ¹⁸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, δ¹³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).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We also thank Dokon Co. Ltd., Taku Ajioka, Ke Wang, Yoshie Nakai, Eriko Noma, and Keiko Ohnishi for help in the sampling, core description, and subsampling, and Stefanie Müller for discussion.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.704332/full#supplementary-material

References

Bronk Ramsey, B. (2008). Deposition Models for Chronological Records. Quat. Sci. Rev. 27 (1–2), 42–60. doi:10.1016/j.quascirev.2007.01.019

CrossRef Full Text | Google Scholar

Bronk Ramsey, C. (2009). Dealing With Outliers and Offsets in Radiocarbon Dating. Radiocarbon 51 (3), 1023–1045. doi:10.1017/s0033822200034093

CrossRef Full Text | Google Scholar

Bronk Ramsey, C. (2021). OxCal 4.4. Available at: http://c14.arch.ox.ac.uk/oxcal (Accessed May 20, 2021).

Google Scholar

Bronk Ramsey, C., and Lee, S. (2013). Recent and Planned Developments of the Program OxCal. Radiocarbon 55 (2), 720–730. doi:10.1017/s0033822200057878

CrossRef Full Text | Google Scholar

Bush, R. T., and McInerney, F. A. (2013). Leaf Wax n-Alkane Distributions in and Across Modern Plants: Implications for Paleoecology and Chemotaxonomy. Geochim. Cosmochim. Acta 117, 161–179. doi:10.1016/j.gca.2013.04.016

CrossRef Full Text | Google Scholar

Chen, X.-Y., Blockley, S. P. E., Staff, R. A., Xu, Y.-G., and Menzies, M. A. (2021). Improved Age Estimates for Holocene Ko-g and Ma-f∼j Tephras in Northern Japan Using Bayesian Statistical Modelling. Quat. Geochronol. [Epub ahead of print].

Google Scholar

Chen, X.-Y., Blockley, S. P. E., Tarasov, P. E., Xu, Y.-G., McLean, D., Tomlinson, E. L., et al. (2016). Clarifying the Distal to Proximal Tephrochronology of the Millennium (B-Tm) Eruption, Changbaishan Volcano, Northeast China. Quat. Geochronol. 33, 61–75. doi:10.1016/j.quageo.2016.02.003

CrossRef Full Text | Google Scholar

Chen, X.-Y., McLean, D., Blockley, S. P. E., Tarasov, P. E., Xu, Y.-G., and Menzies, M. A. (2019). Developing a Holocene Tephrostratigraphy for Northern Japan Using the Sedimentary Record From Lake Kushu, Rebun Island. Quat. Sci. Rev. 215, 272–292. doi:10.1016/j.quascirev.2019.05.017

CrossRef Full Text | Google Scholar

Chikaraishi, Y., Naraoka, H., and Poulson, S. R. (2004). Hydrogen and Carbon Isotopic Fractionations of Lipid Biosynthesis Among Terrestrial (C3, C4 and CAM) and Aquatic Plants. Phytochemistry 65, 1369–1381. doi:10.1016/j.phytochem.2004.03.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Dansgaard, W. (1964). Stable Isotopes in Precipitation. Tellus 16, 436–468. doi:10.3402/tellusa.v16i4.8993

CrossRef Full Text | Google Scholar

De Jonge, C., Hopmans, E. C., Zell, C. I., Kim, J.-H., Schouten, S., and Sinninghe Damsté, J. S. (2014). Occurrence and Abundance of 6-Methyl Branched Glycerol Dialkyl Glycerol Tetraethers in Soils: Implications for Palaeoclimate Reconstruction. Geochim. Cosmochim. Acta 141, 97–112. doi:10.1016/j.gca.2014.06.013

CrossRef Full Text | Google Scholar

Eglinton, G., and Hamilton, R. J. (1967). Leaf Epicuticular Waxes. Science 156, 1322–1335. doi:10.1126/science.156.3780.1322

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, J., Sharkey, T., Berry, J., and Farquhar, G. (1986). Carbon Isotope Discrimination Measured Concurrently With Gas Exchange to Investigate CO2 Diffusion in Leaves of Higher Plants. Funct. Plant Biol. 13, 281–292. doi:10.1071/pp9860281

CrossRef Full Text | Google Scholar

Farquhar, G., and Richards, R. (1984). Isotopic Composition of Plant Carbon Correlates With Water-Use Efficiency of Wheat Genotypes. Funct. Plant Biol. 11, 539–552. doi:10.1071/pp9840539

CrossRef Full Text | Google Scholar

Gao, L., Guimond, J., Thomas, E., and Huang, Y. (2015). Major Trends in Leaf Wax Abundance, δ2H and δ13C Values Along Leaf Venation in Five Species of C3 Plants: Physiological and Geochemical Implications. Org. Geochem. 78, 144–152. doi:10.1016/j.orggeochem.2014.11.005

CrossRef Full Text | Google Scholar

Hoffmann, B., Kahmen, A., Cernusak, L. A., Arndt, S. K., and Sachse, D. (2013). Abundance and Distribution of Leaf Wax n-Alkanes in Leaves of Acacia and Eucalyptus Trees Along a Strong Humidity Gradient in Northern Australia. Org. Geochem. 62, 62–67. doi:10.1016/j.orggeochem.2013.07.003

CrossRef Full Text | Google Scholar

Hopmans, E. C., Schouten, S., Pancost, R. D., van der Meer, M. T. J., and Sinninghe Damsté, J. S. (2000). Analysis of Intact Tetraether Lipids in Archaeal Cell Material and Sediments by High Performance Liquid Chromatography/atmospheric Pressure Chemical Ionization Mass Spectrometry. Rapid Commun. Mass Spectrom. 14, 585–589. doi:10.1002/(sici)1097-0231(20000415)14:7<585::aid-rcm913>3.0.co;2-n

PubMed Abstract | CrossRef Full Text | Google Scholar

Huguet, C., Hopmans, E. C., Febo-Ayala, W., Thompson, D. H., Sinninghe Damsté, J. S., and Schouten, S. (2006). An Improved Method to Determine the Absolute Abundance of Glycerol Dibiphytanyl Glycerol Tetraether Lipids. Org. Geochem. 37, 1036–1041. doi:10.1016/j.orggeochem.2006.05.008

CrossRef Full Text | Google Scholar

Igarashi, Y. (2013). Holocene Vegetation and Climate on Hokkaido Island, Northern Japan. Quat. Int. 290‐291, 139–150.

Inagaki, M., Yamamoto, M., Igarashi, Y., and Ikehara, K. (2009). Biomarker Records From Core GH02-1030 Off Tokachi in the Northwestern Pacific Over the Last 23,000 Years: Environmental Changes during the Last Deglaciation. J. Oceanogr. 65, 847–858. doi:10.1007/s10872-009-0070-4

CrossRef Full Text | Google Scholar

Isono, D., Yamamoto, M., Irino, T., Oba, T., Murayama, M., Nakamura, T., et al. (2009). The 1500-year Climate Oscillation in the Midlatitude North Pacific During the Holocene. Geology 37, 591–594. doi:10.1130/g25667a.1

CrossRef Full Text | Google Scholar

Jouzel, J., Koster, R. D., Suozzo, R. J., and Russell, G. L. (1994). Stable Water Isotope Behavior during the Last Glacial Maximum: A General Circulation Model Analysis. J. Geophys. Res. 99 (D12), 25791–25801. doi:10.1029/94jd01819

CrossRef Full Text | Google Scholar

Kawamura, R., Sugi, M., Kayahara, T., and Sato, N. (1998). Recent Extraordinary Cool and Hot Summers in East Asia Simulated by an Ensemble Climate Experiment. J. Meteorol. Soc. Japan 76, 597–617. doi:10.2151/jmsj1965.76.4_597

CrossRef Full Text | Google Scholar

Keeling, R. F., Graven, H. D., Welp, L. R., Resplandy, L., Bi, J., Piper, S. C., et al. (2017). Atmospheric Evidence for a Global Secular Increase in Carbon Isotopic Discrimination of Land Photosynthesis. Proc. Natl. Acad. Sci. U.S.A. 114, 10361–10366. doi:10.1073/pnas.1619240114

PubMed Abstract | CrossRef Full Text | Google Scholar

Koizumi, I., Tada, R., Narita, H., Irino, T., Aramaki, T., Oba, T., et al. (2006). Paleoceanographic History Around the Tsugaru Strait Between the Japan Sea and the Northwest Pacific Ocean since 30 cal kyr BP. Palaeogeogr. Palaeoclimatol. Palaeoecol. 232, 36–52. doi:10.1016/j.palaeo.2005.09.003

CrossRef Full Text | Google Scholar

Leipe, C., Müller, S., Hille, K., Kato, H., Kobe, F., Schmidt, M., et al. (2018). Vegetation Change and Human Impacts on Rebun Island (Northwest Pacific) over the Last 6000 Years. Quat. Sci. Rev. 193, 129–144. doi:10.1016/j.quascirev.2018.06.011

CrossRef Full Text | Google Scholar

Li, X., Sugimoto, A., and Ueta, A. (2017). Spatial and Temporal Variations of Stable Isotopes in Precipitation in Midlatitude Coastal Regions. Hydrol. Process. 31, 3029–3044. doi:10.1002/hyp.11222

CrossRef Full Text | Google Scholar

Minoshima, K., Kawahata, H., and Ikehara, K. (2007). Changes in Biological Production in the Mixed Water Region (MWR) of the Northwestern North Pacific during the Last 27 Kyr. Palaeogeogr. Palaeoclimatol. Palaeoecol. 254, 430–447. doi:10.1016/j.palaeo.2007.06.022

CrossRef Full Text | Google Scholar

Müller, S., Schmidt, M., Kossler, A., Leipe, C., Irino, T., Yamamoto, M., et al. (2016). Palaeobotanical Records From Rebun Island and Their Potential for Improving the Chronological Control and Understanding Human-Environment Interactions in the Hokkaido Region, Japan. The Holocene 26, 1646–1660. doi:10.1177/0959683616641738

CrossRef Full Text | Google Scholar

Nitta, T. (1987). Convective Activities in the Tropical Western Pacific and Their Impact on the Northern Hemisphere Summer Circulation. J. Meteorol. Soc. Japan 65, 373–390. doi:10.2151/jmsj1965.65.3_373

CrossRef Full Text | Google Scholar

Oppenheimer, C., Wacker, L., Xu, J., Galván, J. D., Stoffel, M., Guillet, S., et al. (2017). Multi-proxy Dating the 'Millennium Eruption' of Changbaishan to Late 946 CE. Quat. Sci. Rev. 158, 164–171. doi:10.1016/j.quascirev.2016.12.024

CrossRef Full Text | Google Scholar

Patwardhan, A. P., and Thompson, D. H. (1999). Efficient Synthesis of 40- and 48-membered Tetraether Macrocyclic Bisphosphocholines. Org. Lett. 1, 241–244. doi:10.1021/ol990567o

PubMed Abstract | CrossRef Full Text | Google Scholar

Reimer, P. J., Austin, W. E. N., Bard, E., Bayliss, A., Blackwell, P. G., Bronk Ramsey, C., et al. (2020). The IntCal20 Northern Hemisphere Radiocarbon Age Calibration Curve 0–55 CAL kBP). Radiocarbon 62, 725–757.

Google Scholar

Sakurai, H., Yamamoto, M., Seki, O., Omori, T., and Sato, T. (2021). Cellulose Oxygen Isotopes of Sphagnum and Vascular Plants in a Peat Core Reveal Climate Change in Northern Japan over the Past 2,000 Years. Geochem. Geophys. Geosyst. [Epub ahead of print].

Google Scholar

Schmidt, M., Leipe, C., Becker, F., Goslar, T., Hoelzmann, P., Mingram, J., et al. (2019). A Multi-Proxy Palaeolimnological Record of the Last 16,600 Years from Coastal Lake Kushu in Northern Japan. Palaeogeogr. Palaeoclimatol. Palaeoecol. 514, 613–626. doi:10.1016/j.palaeo.2018.11.010

CrossRef Full Text | Google Scholar

Schouten, S., Hopmans, E. C., Rosell-Melé, A., Pearson, A., Adam, P., Bauersachs, T., et al. (2013). An Interlaboratory Study of TEX86and BIT Analysis of Sediments, Extracts, and Standard Mixtures. Geochem. Geophys. Geosyst. 14, 5263–5285. doi:10.1002/2013gc004904

CrossRef Full Text | Google Scholar

Schouten, S., Hopmans, E. C., and Sinninghe Damsté, J. S. (2012). The Organic Geochemistry of Glycerol Dialkyl Glycerol Tetraether Lipids: A Review. Org. Geochem. 54, 19–61. doi:10.1016/j.orggeochem.2012.09.006

CrossRef Full Text | Google Scholar

Schubert, B. A., and Jahren, A. H. (2015). Global Increase in Plant Carbon Isotope Fractionation Following the Last Glacial Maximum Caused by Increase in Atmospheric pCO2. Geology 43, 435–438. doi:10.1130/g36467.1

CrossRef Full Text | Google Scholar

Sinninghe Damsté, J. S. (2016). Spatial Heterogeneity of Sources of Branched Tetraethers in Shelf Systems: The Geochemistry of Tetraethers in the Berau River delta (Kalimantan, Indonesia). Geochim. Cosmochim. Acta 186, 13–31. doi:10.1016/j.gca.2016.04.033

CrossRef Full Text | Google Scholar

Sugimoto, A., Higuchi, K., and Kusakabe, M. (1988). Relationship between δD and δ18O Values of Falling Snow Particles from a Separate Cloud. Tellus 40B, 205–213. doi:10.1111/j.1600-0889.1988.tb00290.x

CrossRef Full Text | Google Scholar

Tipple, B. J., and Pagani, M. (2013). Environmental Control on Eastern Broadleaf forest Species' Leaf Wax Distributions and D/H Ratios. Geochim. Cosmochim. Acta 111, 64–77. doi:10.1016/j.gca.2012.10.042

CrossRef Full Text | Google Scholar

Yamamoto, M., Oba, T., Shimamune, J., and Ueshima, T. (2004). Orbital-Scale Anti-Phase Variation of Sea Surface Temperature in Mid-Latitude North Pacific Margins During the Last 145,000 Years. Geophys. Res. Lett. 31, L16311. doi:10.1029/2004gl020138

CrossRef Full Text | Google Scholar

Yamamoto, M. (2009). Response of Mid-Latitude North Pacific Surface Temperatures to Orbital Forcing and Linkage to the East Asian Summer Monsoon and Tropical Ocean-Atmosphere Interactions. J. Quat. Sci. 24, 836–847. doi:10.1002/jqs.1255

CrossRef Full Text | Google Scholar

Yamamoto, M., Suemune, R., and Oba, T. (2005). Equatorward Shift of the Subarctic Boundary in the Northwestern Pacific During the Last Deglaciation. Geophys. Res. Lett. 32, L05609. doi:10.1029/2004gl021903

CrossRef Full Text | Google Scholar