The Japan Sea was a semi-closed marginal sea mainly connected to the subarctic northwestern Pacific via shallow seaways during the late Miocene. We use a multiple regression analysis with common extant radiolarian species groups to estimate the sea-surface temperature (SST) for the period between 9.1 and 5.3 Ma. Our results show a cooling of 8 °C between 7.9 and 6.6 Ma, when the SST dropped from 24 °C to 16 °C. We infer that this cooling dominantly reflects wintertime cooling related to an intensified East Asian winter monsoon. On the other hand, cooling of the summertime SST occurred from 6.6 to 5.8 Ma, suggesting that the late Miocene global cooling is composed of a wintertime cooling phase from 7.9 to 6.6 Ma and summertime cooling phase from 6.6 to 5.8 Ma.

During the late Miocene, fundamental changes in the Earth’s terrestrial ecosystems occurred between 8 and 5 Ma, including an expansion of C4 grassland and the aridification of Asia (e.g., Strömberg, 2011). The reconstruction of sea-surface temperature (SST) variations, based on the organic Uk37 proxy (defined by the relative abundances of C37 alkenones), revealed a sustained cooling event at middle to high latitudes in both hemispheres (Ocean Drilling Program [ODP] Sites 1028, 883, 884, and 887) (Herbert et al., 2016). This cooling event, which is known as the late Miocene global cooling (LMGC), spans the period from 7.8 to 5.8 Ma (Herbert et al., 2016). This interval was also characterized by a sustained worldwide decrease in the inorganic δ13C signature of oceanic deep-water masses (e.g., Holbourn et al., 2018). Holbourn et al. (2018) suggested an intensification of the East Asian winter monsoon (EAWM), which probably intensified marine primary productivity and thus resulted in a higher accumulation of organic carbon in the deep sea because of a steepening in the surface-to-deep organic δ13C gradient, as recorded in the South China Sea (ODP Site 1146). However, the evolution of the EAWM and the impact of the LMGC on northeastern Asia and the northwestern Pacific remain poorly constrained.

The EAWM is driven by the thermal contrast during winter between the Siberian High, which develops over the northern Asian continent, and the Aleutian Low, which forms over the high latitudes of the North Pacific; as a result, cold, dry northwesterly winds blow from Siberia across the Japan Sea (Tada et al., 2016). Thus, the Japan Sea is a strategic location for monitoring the dynamics of the EAWM over northeastern Asia. The Japan Sea is presently a marginal sea with limited oceanic seawater exchange with the East China Sea and the North Pacific through channels that are very shallow (<130 m) and narrow (Fig. 1) (Gamo et al., 2014). In contrast, despite a progressive uplift of northeastern Japan between 9 and 3 Ma, documented by field geological research (e.g., Nakajima et al., 2006; Kharakhinov, 2010), the Japan Sea was open to the North Pacific during the late Miocene, and exchanges of deep water could have occurred between the Japan Sea and the North Pacific (Kozaka et al., 2018). Reconstruction of SST in the Japan Sea during the late Miocene can thus provide valuable information about the long-term evolution of the EAWM. SST variations in the northwestern Pacific are closely related to EAWM dynamics because a strong EAWM strengthens the Pacific Walker circulation and likely weakens the El Niño–Southern Oscillation (e.g., Wang et al., 2005). Calcareous microfossils that are useful for SST estimation are rare in the Japan Sea because of its shallow calcite-compensation depth (e.g., Lee et al., 2000). Instead, the late Miocene sediments in the Japan Sea contain abundant well-preserved diatoms and radiolarians, which have been used for paleoceanographic reconstructions (Tada, 1994; Matsuzaki et al., 2018). Radiolarians are holoplanktonic Protista that are sensitive to changes in SST (Suzuki and Not, 2015), and they have been used to reconstruct Pleistocene SST worldwide with an error of ∼1.0 °C (Hernández-Almeida et al., 2017; Matsuzaki et al., 2019).

In this study, we propose a SST reconstruction method for the Miocene-Pliocene that uses extant radiolarian species whose biogeographic provinces at that time were comparable to those of today, applying a multiple regression analysis to constrain the evolution of the EAWM in the Japan Sea between 8 and 5 Ma.

We analyzed radiolarian assemblages in 116 samples collected by International Ocean Discovery Program (IODP) Expedition 346 from Site U1425, which is on the Yamato Bank (39°29.44′N, 134° 26.55′E) at a water depth of 1909 m (Fig. 1). We used the composite depth scale established by Irino et al. (2018) (CCSF-D_patched version 2). We followed the protocol detailed in Matsuzaki et al. (2018) for the treatment of sediment. The data of these samples were combined with the data of 71 samples from the same site analyzed by Matsuzaki et al. (2018), providing us a data set of 187 samples spanning the period from 9.1 to 5.3 Ma (see the Supplemental Material1). We used the age model of Kurokawa et al. (2019), who tuned the gamma-ray attenuation records from Site U1425 to the short-eccentricity cycle (100 yr) of Laskar et al. (2004), assuming no phase lags. According to this age model, the sampling interval of the collected samples corresponds to ∼10–30 k.y.

We selected extant radiolarian species whose Miocene biogeographic provinces were close to those of today by using the modern data set of Matsuzaki and Itaki (2017), which covers both the low and high latitudes of the northwestern Pacific. The method used is detailed in the Supplemental Material. Before conducting the regression analysis, we verified that the relative abundances of the selected species have Gaussian distributions with annual SST obtained from the World Ocean Atlas 2013 (https://www.nodc.noaa.gov/OC5/woa13/woa13data.html; Locarnini et al., 2013) in the northwestern Pacific. Thus, we performed a normal Q-Q plot, which is a graphical method of comparing two probability distributions by plotting their quantiles against each other. For each species that showed a non-Gaussian distribution, a log10 transformation or a Box-Cox transformation with λ = –0.5 (square root normalization) was performed. Species showing a non-Gaussian distribution after these two transformations were discarded. Then, we checked multicollinearity among the remaining species to identify intercorrelations among them, and grouped together those species showing a high intercorrelation (Jöreskog et al., 2016). Lastly, we conducted a sophisticated multiple regression analysis using R (version 3.5.2). To ensure meaningful multiple regression analysis results, we performed diagnostic linear regression plots. We obtained the following SST equation with an adjusted R2 of 0.84 and a standard error of 2.1 °C. The individual species are represented in relative abundances (%).

The reconstructed median (quartile Q2) annual SST between 9.1 and 5.3 Ma at IODP Site U1425 is 18.9 °C (N = 187) (Fig. 2C). The first- and third-quartile (Q1 and Q3, respectively) SSTs are 16.9 °C and 20.5 °C, respectively. The interquartile range (IQR) is ∼3.6 °C. At ODP Site 1208, at almost the same latitude as IODP Site U1425 in the northwestern Pacific (Fig. 1), the median alkenone-based SST is 23.1 °C during the same time interval (N = 84) (Fig. 2D). The IQR for that interval at ODP Site 1208 is 2.0 °C (Q1 = 21.9 °C and Q3 = 23.9 °C). At the high-latitude ODP Sites 883 and 884 in the northwestern Pacific (Fig. 1), the median alkenone-based SST is 9.5 °C (N = 277) and the IQR is 3.5 °C (Q1 = 8.4 °C and Q3 = 11.9 °C) (Fig. 2E). Thus, the median SST in the Japan Sea (18.9 °C) during the studied interval is 4.2 °C lower than the median recorded at ODP Site 1208 and 9.4 °C higher than that recorded at ODP Sites 883 and 884. The IQRs at IODP Site U1425 and at ODP Sites 883 and 884, 3.6 °C and 3.5 °C, respectively, are similar, whereas the estimated IQR at ODP Site 1208, 2 °C, is lower. Indeed, between 7.9 and 6.6 Ma, the annual SST at IODP Site U1425 and at ODP Sites 883 and 884 decreased similarly, from ∼24.0 °C to 16.0 °C, and from 16.0 °C to 8.0 °C, respectively. In contrast, at ODP Site 1208, SST decreased by 4 °C between 6.6 and 5.8 Ma.

Our data indicate that the Japan Sea experienced sustained cooling between 7.9 and 6.6 Ma, when the annual SST dropped from 24.0 °C to 16.0 °C (Fig. 2C). Synchronous cooling of similar magnitude (8.0 °C) is recorded in the subarctic North Pacific at ODP Sites 883 and 884, where the SST decreased from 16.0 °C to 8.0 °C between 7.9 and 6.6 Ma (Herbert et al., 2016) (Figs. 2C and 2E). This similarity is possibly caused by closure of the eastern seaway in response to the tectonic uplift of Honshu island (Japan) between 9 and 3 Ma (Tada, 1994), which caused the Japan Sea to be connected only to the subarctic northwestern Pacific starting at ca. 7–8 Ma (Kharakhinov, 2010; Matsuzaki et al., 2018). Therefore, the Japan Sea was probably affected more by subarctic waters flowing southward from the Sea of Okhotsk during the late Miocene, whereas ODP Site 1208 stayed to the south of the subpolar front. In addition, in the modern mid- to high-latitude northwestern Pacific, the radiolarian flux is higher from autumn to winter, with the maximum flux >1200 skeletons m−2 d−1 in December, whereas during the summer, the radiolarian flux is ∼200 skeletons m−2 d−1 (Itaki et al., 2008). Therefore, it is reasonable to infer that the radiolarian-based annual mean SST is highly influenced by the winter SST. On the other hand, alkenone-based SSTs are influenced by bloom season of the haptophytes, when the highest alkenone fluxes are recorded (e.g., Jonas et al., 2017). At high latitudes of the North Pacific, the bloom season is likely from autumn to winter (e.g., Harada et al., 2003; Seki et al., 2007). Thus, similar to the radiolarian-based SST, the alkenone-based SST in the subarctic North Pacific likely reflects autumn–winter conditions. Therefore, the synchronicity of the cooling between IODP Site U1425 and ODP Sites 883 and 884 and the similarity in their cooling magnitude suggest that wintertime SST in the northwestern Pacific decreased substantially from 7.9 to 6.6 Ma. Herbert et al. (2016) proposed that a LMGC event occurred between 7.8 and 5.8 Ma at middle to high latitudes of both hemispheres. Therefore, the winter cooling event recorded in the Japan Sea and the subarctic northwestern Pacific between 7.9 and 6.6 Ma is likely related to the LMGC, and the local paleogeographic configuration in which the Japan Sea was isolated from the mid-latitude northwestern Pacific Ocean amplified the seasonal cooling of the Japan Sea during winter.

Although a decrease in SST is recorded worldwide during the LMGC, the deep-sea marine δ18O signature does not imply a corresponding increase in the volume of the Earth’s major ice sheets between 7.9 and 5.8 Ma (Fig. 2A) (Holbourn et al., 2018). Moreover, when SSTs of middle to high latitudes of both hemispheres are compared, the magnitude of the SST cooling during the LMGC seems to be larger at high latitudes of the Northern Hemisphere, probably because of larger land area in the Northern Hemisphere. Cooling of the high latitude Northern Hemisphere intensified the Siberian High, which in turn enhanced EAWM and further cooling of the northwestern Pacific. Clay mineral and palynological studies indicated that the climate in East Asia became colder and dryer as a result of the probable intensification of the EAWM between 8.0 and ca. 6.5 Ma, possibly related to the uplifting of the Tibetan Plateau (e.g., An et al., 2001; Wan et al., 2007; Miao et al., 2017; Clift et al., 2019). The inferred low wintertime SST in the northwestern Pacific between 7.9 and 6.6 Ma appears thus to be associated with the EAWM intensification.

We used the method proposed in this study to reconstruct radiolarian-based wintertime SST with the data set of Matsuzaki et al. (2019) from the East China Sea (ECS), an area influenced by the EAWM (Tada et al., 2016), and compared the results with mean grain size of quartz data from the Chinese Loess Plateau (Sun et al., 2006), which is an EAWM proxy (Sun et al., 2006). The result shows that during the last 150 k.y., the radiolarian-based SST and the grain-size data have varied in a similar manner, with the larger quartz grain-size (i.e., strong winter monsoon) periods corresponding to the periods of low wintertime SST (Fig. 3). Therefore, the variability of wintertime SST observed at IODP Site U1425 is likely related to the EAWM. Although there are no published quartz grain-size data for the interval between 8 and 7 Ma, the quartz grain size in the Chinese Loess Plateau seems to be large and sensitive to the obliquity cycle between 7 and 6 Ma (Sun et al., 2010). Thus, we attribute the drastic decrease in wintertime SST between 7.9 and 6.6 Ma to a steady intensification of the EAWM affecting the climate throughout East Asia during the LMGC, probably because of an intensification of the Siberian High.

The decrease in wintertime SST at IODP Site U1425 and ODP Sites 883 and 884 is synchronous with the sharp negative shift in the deep-sea δ13C at ODP Site 1146 (Fig. 2B) and also with the expansion of C4 vegetation (Strömberg, 2011; Herbert et al., 2016; Holbourn et al., 2018). The C4 vegetation has low water losses and sequesters more organic carbon during its photosynthesis (Edwards and Walker, 1983). Thus, although speculative, the expansion of C4 vegetation from 7.9 to 6.6 Ma and the consequent enhancement of organic carbon sequestration and burial could have resulted in pCO2 decrease and global cooling and aridification.

The SST changes at ODP Site 1208 show a pattern different from those at IODP Site U1425 and ODP Sites 883 and 884 (Figs. 2C–2E). At ODP Site 1208, which is located at similar latitude to IODP Site U1425, a cooling of 1 °C is recorded between 7.9 and 6.6 Ma (Fig. 2D), whereas a cooling of 4 °C is observed between 6.6 and 5.8 Ma during the latter half of the LMGC (Fig. 2D). In middle latitudes of the northwestern Pacific, the bloom season of the haptophytes that produce alkenones is from late spring to summer, with fluxes of as much as 16.5 μg m−2 d−1 (Sawada et al., 1998). Thus, in the mid-latitudes of the northwestern Pacific, alkenone-based SST is highly influenced by summer SST. If this interpretation is correct, our data suggest that the earlier half of the LMGC during 7.9–6.6 Ma was triggered by a steady decrease in the wintertime SST of 8 °C from 7.9 to 6.6 Ma in association with an intensification of EAWM and changes in on-land ecosystems and carbon cycle. Summertime SST seems to be less affected by the intensification of the EAWM from 7.9 to 6.6 Ma, with an SST decrease of only 1 °C, whereas summer SST cooling occurred during 6.6–5.8 Ma with slight wintertime SST warming of ∼2–3 °C. Therefore, the time period between 7.9 and 6.6 Ma was likely characterized by distinct winter cooling with high seasonality, probably due to intensification of EAWM. On the other hand, summertime SST recorded a 4 °C cooling from 6.6 Ma to 5.8 Ma, while wintertime SST probably slightly increased, which resulted in a decrease in SST seasonality.

We reconstructed, for the first time, annual SST in the Japan Sea during the late Miocene with an error of 2.1 °C (R2 = 0.84) by using data of radiolarian species at IODP Site U1425. The reconstructed SSTs are inferred to be wintertime SSTs; thus, our reconstruction indicates that the Japan Sea experienced wintertime cooling of 8 °C between 7.9 and 6.6 Ma because of an intensified EAWM in East Asia. The summertime SSTs were probably were less affected, and an SST decrease of 1 °C is observed from 7.9 and 6.6 Ma, which resulted in high seasonality in the SST. On the other hand, summertime SST of the northwestern Pacific started to cool ca. 6.6 Ma, causing a decreasing in the seasonality of the SST.

We thank Shiming Wan and two anonymous reviewers for their helpful comments. This work was supported by Japan Society for the Promotion of Science grant 19KK0089 to K. Matsuzaki, grants JPMXS05R2900001 and 19H05595 to M. Yamamoto, and grant 18H01279 to T. Sagawa. We thank the IODP core curators at the Kochi Core Center (Nankoku, Japan) for providing samples. We thank A. Holbourn, W. Kuhnt, M. Ikeda, and T. Itaki for advice.

1Supplemental Material. Details of the method used to estimate radiolarian-based sea-surface temperature. Please visit https://doi.org/10.1130/GEOL.S.12331130 to access the supplemental material, and contact [email protected] with any questions.