Articles | Volume 16, issue 1
Research article
 | Highlight paper
07 Feb 2020
Research article | Highlight paper |  | 07 Feb 2020

Modal shift in North Atlantic seasonality during the last deglaciation

Geert-Jan A. Brummer, Brett Metcalfe, Wouter Feldmeijer, Maarten A. Prins, Jasmijn van 't Hoff, and Gerald M. Ganssen

Changeover from a glacial to an interglacial climate is considered as transitional between two stable modes. Palaeoceanographic reconstructions using the polar foraminifera Neogloboquadrina pachyderma highlight the retreat of the Polar Front during the last deglaciation in terms of both its decreasing abundance and stable oxygen isotope values (δ18O) in sediment cores. While conventional isotope analysis of pooled N. pachyderma and G. bulloides shells shows a warming trend concurrent with the retreating ice, new single-shell measurements reveal that this trend is composed of two isotopically different populations that are morphologically indistinguishable. Using modern time series as analogues for interpreting downcore data, glacial productivity in the mid-North Atlantic appears limited to a single maximum in late summer, followed by the melting of drifting icebergs and winter sea ice. Despite collapsing ice sheets and global warming during the deglaciation, a second “warm” population of N. pachyderma appears in a bimodal seasonal succession, separated by the subpolar G. bulloides. This represents a shift in the timing of the main plankton bloom from late to early summer in a “deglacial” intermediate mode that persisted from the glacial maximum until the start of the Holocene. When seawater temperatures exceeded the threshold values, first the “cold” (glacial) then the “warm” (deglacial) populations of N. pachyderma disappeared, whilst G. bulloides with a greater tolerance to higher temperatures persisted throughout the Holocene to the present day in the midlatitude North Atlantic. Single-specimen δ18O of polar N. pachyderma reveals a steeper rate of ocean warming during the last deglaciation than appears from conventional pooled δ18O average values.

1 Introduction

1.1 Seasonality and single-foraminiferal analysis (SFA)

Stable oxygen isotopes (δ18O) of pooled foraminifera have been used as key tracers of water masses (e.g. Epstein and Mayeda, 1953; Frew et al., 2000), ice-volume and sea-level fluctuations (e.g. Grant et al., 2012; Waelbroeck et al., 2002; Shackleton, 1987) over glacial–interglacial cycles (e.g. Pearson, 2012; Waelbroeck et al., 2005). Technical advances allow (Killingley et al., 1981; Schiffelbein and Hills, 1984; Oba, 1990, 1991; Billups and Spero, 1996) for the routine analysis of the stable isotopic composition of single microscopic shells (Feldmeijer et al., 2015; Lougheed et al., 2018; Metcalfe et al., 2015, 2019; Frizt-Endres et al., 2019; van Sebille et al., 2015) and chambers (Takagi et al., 2015, 2016; Lougheed et al., 2018; Pracht et al., 2019) of foraminifera permitting the resolution of (sub)seasonal contrasts in seawater temperature (Ganssen et al., 2011; Wit et al., 2010; Groenveld et al., 2019), in lieu of pooled specimens that capture an averaged state of the system on longer timescales. The use of single-shell oxygen isotope analysis allows for moving beyond the “average” state of the climate system, as expressed in pooled specimen analysis to observe the interspecimen variance (e.g. Leduc et al., 2009; Koutavas et al., 2006; Koutavas and Joanides, 2012; Scussolini et al., 2013) that includes seasonal differences (Feldmeijer et al., 2015; Ganssen et al., 2011; Metcalfe et al., 2015; Wit et al., 2010). Seasonal changes during these glacial–interglacial cycles have rarely been addressed, although resolving seasonal contrasts would significantly improve our understanding of past climate change (Huybers, 2006; Schmittner et al., 2011).

1.2 Aims and objectives

As the largest ocean carbon sink in the Northern Hemisphere, the North Atlantic Ocean (Gruber et al., 2002) exhibits strongly seasonal productivity in the present day. Deep wind-driven mixing in winter resupplies the photic zone with nutrients brought up from subsurface depths (Falkowski and Oliver, 2007), leading to phytoplankton blooms and maxima in the abundance of zooplankton including planktonic foraminifera during the onset of summer stratification, followed by a decrease as oligotrophic summer conditions develop. Present-day temperature conditions in the midlatitude North Atlantic (Fig. 1) preclude the occurrence of Neogloboquadrina pachyderma (Kretschmer et al., 2016), the species being restricted to the (sub)polar water masses in the high-latitude North Atlantic (Kohfeld et al., 1996). With the southward shift of the Polar Front that accompanied the last glacial, more favourable conditions developed, whereas Globigerina bulloides (Ganssen and Kroon, 2000) existed throughout. Here, we analyse single-shell stable oxygen and carbon isotopes of the planktonic foraminifera left-coiling N. pachyderma (e.g. Mikis et al., 2019; Metcalfe et al., 2019) and G. bulloides (e.g. Ganssen et al., 2011; Metcalfe et al., 2015) in a sediment core from the Iceland Basin in the midlatitude Atlantic in order to address the direction of mean annual temperature change and seasonal changes during the past deglaciation. Given that present conditions in the mid-North Atlantic are an anathema to the polar species N. pachyderma, this species is only used during the glacial and deglacial sections of the core.

Figure 1Map of North Atlantic study area with location of core T88-3P and the position of the sediment trap time series from the (A) Iceland Plateau (IP), (B) Irminger Sea (IRM) and (C) North Atlantic Bloom Experiment (NABE). Base maps represent the sea surface temperature for January–February–March of the Last Glacial Maximum, based on the Multiproxy Approach for the Reconstruction of the Glacial Ocean Surface (MARGO) database (Kučera et al., 2005), and the modern ocean, based on World Ocean Atlas 1998 (as used by MARGO).

2 Methodology

Piston core T88-3P (56.49 N, 27.80 W; Fig. 1) was taken on the eastern flank of the Mid-Atlantic Ridge during the 1988 RV Tyro expedition of the Actuomicropalaeontology Palaeoceanography North Atlantic Project (APNAP) II. Piston core T88-3P measures 937 cm in length (Figs. 2 to 4) and was retrieved from above both the modern and glacial carbonate compensation depth (CCD) (core water depth: 2819 m), ensuring minimal bias by carbonate dissolution. Core sections were manually split into a working half and an archive half.

2.1 X-ray fluorescence core scanning and composite images

Archive halves of each section of the entire piston core were analysed at 1 cm downcore resolution using the Avaatech X-ray fluorescence (XRF) core scanner (Richter et al., 2006) at the Royal NIOZ (Fig. 4). Optical line scanning was first performed on the split halves allowing a detailed and accurate description of visual and chromatic changes in core texture (Supplement Fig. S1). Prior to XRF analysis, the surface of the archive halves was scraped cleaned and each section was covered in SPEXCerti Ultralene® ultra-thin (4 µm) film. Bulk chemical composition was measured using energy-dispersive fluorescence radiation as elemental intensities in counts per second (CPS) at 10 kV (for 10 s) and at 50 kV (for 40 s). Despite limitations on the accuracy and precision (Weltje and Tjallingii, 2008) by matrix effects, sediment (e.g. water content; grain size) and measurement properties (e.g. surface irregularities) as well as machines settings used (outlined above), the reliability for the elements Ca and Ti used herein is well established (Weltje and Tjallingii, 2008). To further minimise error, counts are expressed as log ratios of two elements (Weltje and Tjallingii, 2008). Herein, log(Ca∕Ti) is used as a proxy for two end-members: marine productivity ([Ca]) and detrital terrestrial material ([Ti]), with minor contribution to [Ca] via detrital carbonate, which directly relates to ice-rafted debris (IRD; Fig. 4).

2.2 Abundance counts

The core sections of the entire working half were sampled every centimetre, resulting in 1 cm sample slices that were each washed over a 63 µm sieve mesh, dried overnight at ∼75C and subsequently size fractionated into 63–150 µm and > 150 µm. For abundance counts of planktonic foraminifera, slices every 4 cm were used; the counts were performed on G. bulloides, N. pachyderma, “other foraminifera” and terrigenous grain in the > 150 µm size fraction. The relative abundance was calculated from a total of both planktonic foraminifera (G. bulloides; N. pachyderma; and “other foraminifera”) and IRD (stained quartz; cloudy quartz; bright quartz; quartz; sandstone, igneous, obsidian glass; rhyolitic glass and “other”) counts, and this does have complications for the relative abundance of foraminiferal species as it is a closed sum. Everything was counted and identified on a minimum of 200 grains after splitting with an OTTO microsplitter. The ratio of N. pachyderma to G. bulloides (Fig. 4) is expressed as

(1) Ratio of NPS = N. pachyderma ( N. pachyderma + G. bulloides ) .

2.3 Single-foraminifera stable isotope geochemistry (δ18O; δ13C)

For single-shell stable isotope analysis, a continuous-flow isotope ratio mass spectrometer was used (Feldmeijer et al., 2015; Metcalfe et al., 2015) based on modifications (Breitenbach and Bernasconi, 2011) to the microvolume setup (Spötl and Vennemann, 2003). Slices for isotope analysis were selected first at 10 cm resolution and then at specific sections downcore every 2 cm. Two sections were analysed: a deglacial section between 300 and 420 cm and a glacial section between 515 and 565 cm. For each slice, up to 20 shells of both left-coiling N. pachyderma and G. bulloides were picked at random from the 250–300 µm size fraction (Figs. 5 to 6). No morphological differences were observed among the picked left-coiling N. pachyderma. Each specimen was placed into a 4.5 mL exetainer vial, and the ambient air was replaced by He and subsequently digested in concentrated H3PO4 (45 C for 160 min). The resultant CO2–He gas mixture was transported to the GasBench II using a He flow through a flushing needle system where water was extracted from the gas using a Nafion tubing. The purified CO2 was analysed in a Thermo Finnigan Delta+ mass spectrometer after separation from other gases in a gas chromatography (GC) column. Isotope values are reported in the standard δ denotation with the ratio of heavy to light isotopes (δ18O) in per mil (‰) versus Vienna Pee Dee belemnite (V-PDB). The reproducibility of an international carbonate standard (IAEA-CO1) analysed is < 0.12 ‰ (1σ) for both δ18O and δ13C, measured within the same run and at similar quantities (i.e. producing similar amplitude on mass 44) to a single foraminifer. Based on the amplitude of mass 44, which correlates with shell weight, foraminifera are estimated to weigh > 10 µg. Following Ganssen et al. (2011), data were screened for anonymous values, leading to outlier-corrected values (red data points in Figs. 5b, c, and 6b, c) for both the deglacial (N. pachyderma: ntotal=414; noutlier corrected=388; noutlier=26; G. bulloides: ntotal=439; noutlier corrected=424; noutlier=15) and glacial sections (N. pachyderma: ntotal=490; noutlier corrected=474; noutlier=16; G. bulloides: ntotal=496; noutlier corrected=474; noutlier=22). Due to the nature of the mixture analysis (see Sect. 2.5), the resultant outlier-corrected populations differ from non-outlier-corrected populations – although most of the uncorrected samples bearing outliers were beyond the models capacity to “unmix” the results.

Figure 2Raw pooled foraminifera stable oxygen isotope data and computed averages. (a) Pooled stable oxygen isotope data of G. bulloides and (b) G. glutinata; for each sample, two groups (white and black dots) were measured and (c) the average was computed for both species (lines in a–c).


Figure 3Age models of core T88-3P and their respective sediment accumulation rates (SARs). (a) Radiocarbon age model (see Table 1 and Table S1 in the Supplement for radiocarbon ages) and (b) the estimated SAR. (c) Oxygen isotope stratigraphy of pooled measurements of G. glutinata and G. bulloides (see Fig. 2) tuned to the North Greenland Ice Core Project (NGRIP) (see Table S2 in the Supplement for tie-point estimates). For comparison, the radiocarbon age model (blue line) is plotted alongside the tuned oxygen isotope age model (black line). IRD downcore is plotted (grey) alongside the age models. (d) The estimated SAR for the tuned oxygen isotope age model. Green panels in panels (b) and (d) highlight depths in core where single-foraminifera stable isotope analysis was performed. See Table 1 and Table S1 in the Supplement for radiocarbon measurements and Table S2 in the Supplement for δ18O tie points.


Figure 4Sediment and abundance data against depth in core. (a) Average oxygen stable isotopes of G. bulloides (blue) and G. glutinata (red) from Fig. 2. (b) Relative abundance of IRD, calculated as the amount of IRD relative to both foraminifera and IRD from approximately 200 particles. Grey area reflects the relative abundance of IRD, whereas the white area reflects the relative abundance of foraminifera. (c) The logarithmic ratio of Ca and Ti (log(Ca∕Ti)) counts per second (CPS) as measured by X-ray fluorescence (XRF) core scanning. (d) The relative abundance of the planktonic foraminifera species G. bulloides (blue) and N. pachyderma (green); the relative abundance is based on the counts of IRD and foraminifera (see panel b). Depths reflect the midpoint of the sample. (e) Ratio of the relative abundance of N. pachyderma and G. bulloides (see Eq. 1).


Figure 5Single-foraminifera stable isotope data: N. pachyderma. (a) Pooled average oxygen stable isotopes of G. bulloides (blue) and G. glutinata (red), from Fig. 2, and ice-rafted debris (grey) for the same samples (Figs. 3 and 4). Green areas represent the deglacial (b, d) and glacial (c, e) samples. (b–c) Raw stable isotope values (grey) of N. pachyderma for each sample and the sample's respective outliers (red). Overlain are the statistical features of the distribution outlined by a box-and-whisker plot (median, upper and lower quartiles). (d–e) Average values and the standard deviation (blue) of the outlier-corrected oxygen isotope data plotted alongside the abundance of N. pachyderma (green; see Fig. 4) and IRD (grey; see Fig. 4).


Figure 6Single-foraminifera stable isotope data: G. bulloides. (a) Pooled average oxygen stable isotopes of G. bulloides (blue) and G. glutinata (red), from Fig. 2, and ice-rafted debris (grey) for the same samples (Figs. 3 and 4). Green areas represent the deglacial (b, d) and glacial (c, e) samples. (b–c) Raw stable isotope values (grey) of G. bulloides for each sample and the sample's respective outliers (red). Overlain are the statistical features of the distribution outlined by a box-and-whisker plot (median, upper and lower quartiles). (d–e) Average values and the standard deviation (blue) of the outlier-corrected oxygen isotope data plotted alongside the abundance of G. bulloides (green; see Fig. 4) and IRD (grey; see Fig. 4).


2.4 Core stratigraphy and age model

2.4.1 Radiocarbon dating

For radiocarbon dating of core T88-3P, approximately 1 mg of pristine specimens of G. bulloides and N. pachyderma were picked from six samples of core T88-3P and analysed by accelerated mass spectrometry (AMS) at the AMS laboratories of Beta Analytic (Table 1; Fig. 3). The open source MatCal (version 2.0) function for Mathworks MatLab® (Lougheed and Obrochta, 2016) was used to calibrate conventional radiocarbon age to a calendar age, using the Marine13 calibration curve (Reimer et al., 2013) and a reservoir age of 400 14C yr, with an error of 200 14C yr, expressed mathematically as ΔR: 0±20014C yr (Reimer et al., 2013). The 95 % confidence limits for the calendar age, in kyr BP, of each sample are given in Table 1.

Table 1Raw and calibrated radiocarbon ages. Conventional radiocarbon age represents the measured radiocarbon age (see Table S1 in the Supplement) corrected for isotopic fraction.

Download Print Version | Download XLSX

2.4.2 Age model construction

Two age models were produced that used two independent methods: first, a radiocarbon-only age model, and second, δ18O a stratigraphy-only age model (Fig. 3). The first radiocarbon age model uses the six radiocarbon dates, using the maximum likely calendar age (in cal yr BP) from Table 1; the age model consisting of independent age markers places the deglacial period between ∼410 and ∼290 cm downcore. A change in sedimentation occurs ∼10 000 years ago at the core site, from a slow glacial sediment accumulation rate (SAR) (∼10 cm kyr−1) to a rapid interglacial SAR (∼30 cm kyr−1); this may reflect the site's location and the position of the Polar Front (Fig. 1). The radiocarbon date at 500 cm is older than > 35 kyr; therefore, the calibrated age can be considered less robust than the younger radiocarbon dates. A second, independent, downcore age model using δ18O stratigraphy was constructed, as the deeper depths (> 500 cm) of the core are beyond the limits of radiocarbon dating. This stratigraphy utilised the cosmopolitan upper ocean dweller G. glutinata and the subpolar–temperate upper ocean dweller G. bulloides; these species were measured for δ18O and δ13C using pooled specimens (two groups of 5–10 specimens) picked from the 250–300 µm size fraction (Fig. 2), which were placed in monospecific groups within a 15 mL exetainer vial. Analyses followed the same procedures as for single-shell analysis but with considerably better reproducibility of international standards (1σ < 0.10 ‰) for the larger sample mass (∼100µg). The average δ18O of G. bulloides and G. glutinata was tuned to composite record of North Greenland Ice Core Project (NGRIP) (Rasmussen et al., 2008) on the GICC05 timescale (subtracting 50 years to allow for a comparison between BP and b2K). Given the differences in resolution between a marine core and an ice core, the average δ18O of G. bulloides and G. glutinata was tuned to a filtered NGRIP signal: here, we use the average of “envelope” that reproduces the magnitude (highest value, lowest value) using a discrete Fourier transform with a Hilbert finite impulse response (FIR) filter of length 100. The two age models agree for the sections of the core between 0 and 30 000 years; however, after 35 kyr, the two age models diverge. However, this may be simply a function of the limitation of the radiocarbon calibration curve > 35 kyr.

Using two independent age models, two estimates of the sedimentation rate have been made (Fig. 3). The background sedimentation rate varies between 2.5 and 40 cm kyr−1, being noticeably slower during glacial periods and much faster during the Holocene (interglacial). During intervals of high IRD input the sedimentation rate noticeably varies. With respect to intervals chosen for single-foraminiferal analysis, the SAR increases over the first interval (green boxes in Fig. 3) during the deglacial interval, yet SAR stays relatively constant for the second interval chosen (Fig. 3). Although given the depth in core, for this second interval, only one estimate of SAR can be made.

2.5 Statistical analysis: end-member modelling of SFA

Marine sediments reflect an averaged record over time, ranging from months to multiple centuries. However, if the individual components have distinct markers such as different δ18O values, then the original distributions can be statistically unmixed into two or more univariate normally distributed populations using an unmixing function (e.g. Hammer et al., 2001; Weltje, 1997; Weltje and Prins, 2003; Wit et al., 2013). Mixture analysis was carried out on outlier-corrected samples (Fig. 7) using the open-source PAST (version 3.10) palaeontological statistics software (Dempster et al., 1977; Hammer et al., 2001). Using the end-member modelling algorithm of Dempster et al. (1977), PAST estimates the mean, standard deviation and proportion of each population (see Hammer et al., 2001, for a discussion of the assumptions of the mixing model). These solutions can be tested by two methods: the log-likelihood value in which a “better” result produces a less negative value and a minimum in Akaike information criterion (AIC) value indicating that the chosen number of groups has a good fit without subsequent overfitting. An additional output of this mixture analysis is to assign each individual to the most probable population.

2.6 Modern sediment trap record and ocean reanalysis

As the modern analogue of our distinct isotopic “end-members”, we used seasonally resolved sediment trap time series representing the modern polar, subpolar and temperate North Atlantic (Fig. 1). Three such sediment trap records are available (Figs. 1 and 8) from (a) the polar Greenland–Norwegian Sea over the Iceland Plateau (IP; Wolfteich, 1994), (b) the subpolar Irminger Sea (IRM; Jonkers et al., 2010, 2013; Jonkers and Kučera, 2015) and (c) the temperate mid-North Atlantic (NABE48 from the North Atlantic Bloom Experiment; Wolfteich, 1994). Ocean ReAnalysis System 4 (ORAS4) (Balmaseda et al., 2013) was used to complete the temperature and salinity profiles associated with each sediment trap time series, both with respect to time and depth (Fig. 8). Ocean reanalysis data were converted from date into sediment trap cup number using a Mathworks MatLab® function: the monthly temperature and salinity data were first interpolated to 1 d resolution, using the interp1 function; the opening and closing dates of successive cups were then found, temperature and/or salinity presented in figures represent the opening and therefore the closing of the previous cup used. Since the IRM time series represent several years, we generate both a time-averaged flux record as well as an average profile for both temperature and salinity. The time-averaged flux is calculated by finding corresponding bimonthly (cup opening interval: 14 d) trap opening and closing days and averaging the resultant flux.

3 Results

3.1 IRD, abundance

The upper ∼290 cm of core T88-3P is Holocene in age as evidenced by near-uniform values of pooled specimen δ18O values, log(Ca∕Ti), IRD and the ratio of NPS (Fig. 4). Between 290 and 410 cm, the deglacial interval, the ratio of NPS approaches 1.0, IRD approaches 20 % and a minimum log(Ca∕Ti) of 0.9. The minimum in log(Ca∕Ti) occurs prior to the increase in IRD though coeval with the increase in the ratio of NPS. Between 425 and 937 cm, the ratio of NPS and percentage of IRD appear to covary, whilst the log(Ca∕Ti) shows an inverse, with a minimum in log(Ca∕Ti) during IRD events.

3.2 Single-shell δ18O: N. pachyderma

Single-shell analysis of N. pachyderma (Fig. 5) was performed on a deglacial interval (300 to 420 cm) and a glacial interval (515–565 cm). Our glacial results show the abundance during this period has two peaks centred at 515 and 560 cm, with a large proportion of the data occurring within an interval of lower abundance. Single-shell δ18O average values and standard deviation lie between 2.5 and 4.5 ‰, with a remarkable consistent spread in the values; the lightest δ18O values occur during lower species abundance. In comparison, our deglacial results show the abundance decreasing from a peak centred at 380 cm. Single-shell values appear to get heavier from 420 to 360 cm, with a reduced spread, before becoming lighter and having a larger spread between 360 and 340 cm. The δ18O average values and standard deviation during this interval range from 5.5 to 2.5 ‰.

3.3 Single-shell δ18O: G. bulloides

Single-shell analysis of G. bulloides was performed (Fig. 6) on a deglacial interval (300 to 420 cm) and a glacial interval (515–565 cm). Our glacial results show the abundance during this period has a single peaks centred at 540 cm, with a large proportion of the data occurring within an interval of high abundance. Single-shell δ18O average values and standard deviation lie between 3.5 and 1.5 ‰, with a remarkably consistent spread in the values. In comparison, our deglacial results show the abundance decreasing from a peak centred at 420 cm, reaching a lower limit between 380 and 360 cm, before rising again. Single-shell values appear to be relatively consistent from 420 to 360 cm, with a reduced spread, before becoming lighter and having a larger spread between 360 and 340 cm. The δ18O average values and standard deviation during the deglacial interval range from 4 to 1.5 ‰.

3.4 Single-shell δ18O standard deviation

Comparison between the glacial and deglacial samples highlights the change in standard deviation between the two time periods: the glacial samples (515–565 cm) for N. pachyderma have standard deviation corrected for outliers (μ=0.44; min = 0.28; max = 0.61; σ=0.09; ngroups=26; nwithin group=474) and uncorrected for outliers (μ=0.55; min = 0.28; max = 0.85; σ=0.15; ngroups=26; nwithin group=490) lower than the deglacial samples (300–420 cm) either uncorrected for outliers (μ=0.70; min = 0.26; max = 1.52; σ=0.31; ngroups=22; nwithin group=414) or corrected for outliers (μ=0.49; min = 0.13; max = 1.11; σ=0.26; ngroups=22; nwithin group=388). The deglacial interval has a larger range in standard deviation that the glacial intervals, with both the smallest and largest spread as represented by the sample standard deviation. Whereas, the deglacial data for G. bulloides uncorrected (μ=0.46; min = 0.24; max = 0.72; σ=0.13; ngroups=23; nwithin group=439) and corrected (μ=0.38; min = 0.19; max = 0.65; σ=0.13; ngroups=23; nwithin group=424) are somewhat similar to the glacial uncorrected (μ=0.69; min = 0.33; max = 1.53; σ=0.28; ngroups=26; nwithin group=496) and corrected (μ=0.48; min = 0.24; max = 0.88; σ=0.16; ngroups=26; nwithin group=474) data.

3.5 Single-shell δ18O populations

Our results show that δ18O values of both G. bulloides and N. pachyderma are predominately unimodally distributed during the last glacial until about 21 ka (Figs. 5 and 6). In the glacial section (515–565 cm), only a few samples appear to have a second population for N. pachyderma (5 out of 26 samples) and G. bulloides (9 out of 26 samples), with their appearance being more sporadic than systematic. This shows remarkable contrast with the deglacial section (300–420 cm). Whilst the distribution of G. bulloides δ18O values remains predominately unimodal throughout the deglacial section (13 out of 23 samples), the δ18O values of N. pachyderma develop striking bimodality (12 out of 22 samples; Figs. 5 and 7). For N. pachyderma, the distributions can be statistically unmixed into two discrete populations (Hammer et al., 2001) in varying numbers of specimens: one high in δ18O persisting from the glacial (population P1) and a second population low in δ18O appearing at the onset of the deglaciation (population P2). The difference in δ18O between populations P1 and P2 amounts to 0.9±0.4 ‰ and persists for (as estimated by the age models; Fig. 3) about 10 kyr, while absolute values gradually decrease by 1.6 ‰ (Figs. 5 and 7). At the end of the last deglaciation, P1 disappears and the δ18O values of N. pachyderma become once more unimodal, now for P2, shortly before disappearing entirely until the present day. Carbon isotope values (δ13C) measured on the same shells of N. pachyderma do not appear to show this bimodal distribution, precluding the possibility that the two populations in δ18O represent a similar season but grew their shells at different depths given the enrichment and depletion with depth in seawater 13C associated with phytoplankton growth and decay. Since the δ18O values of N. pachyderma exhibit bimodality, our findings downcore could equate with seasonal gradients and species successions, as observed in modern time series from sediment traps deployed in the modern North Atlantic (Figs. 1 and 8) at 48 N (temperate), 59 N (subpolar) and 68 N (polar).

4 Discussion

4.1 Modern analogue

Modern conditions that mimic glacial times downcore are presently found in the polar Greenland–Norwegian Sea where productivity is limited by low light conditions, deep mixing and intermittent sea ice cover (Kučera et al., 2005). At 68 N, late summer insolation and thermal stratification spur a plankton bloom (August–September). At the same time, planktonic foraminifera produce a single high maximum in the shell flux of N. pachyderma with few G. bulloides (Jonkers and Kučera, 2015) at temperatures of 3–5 C, before the arrival of meltwater (Fig. 8a). Further south, at 59 N, in the subpolar Irminger Sea, the flux of N. pachyderma is bimodal, with an early “cold” population being produced in April–May (4–6 C) and a late “warm” population occurring in August–September (7–9 C) that are separated by a single pulse in G. bulloides (Jonkers et al., 2010, 2013) (Fig. 8b). Neither of the N. pachyderma populations from IRM displays significant morphological differences (Jonkers et al., 2010, 2013). By contrast, modern shell fluxes in the temperate North Atlantic at 48 N, close to our core site, are dominated by G. bulloides in early summer, yet completely devoid of N. pachyderma year around (Wolfteich, 1994) (Fig. 8c).

Spatial differences in modern seasonality observed in the polar-to-temperate North Atlantic provide modern analogues for interpreting temporal changes in the sediment record in terms of the seasonal modes developing since the last glacial. During peak glacial times, the northern North Atlantic is covered by sea ice down to 45 N (Kučera et al., 2005) (Fig. 1), except for a short interval in late summer allowing for a period of high productivity dominated by N. pachyderma (P1), as seen in the modern Norwegian–Greenland Sea at 68 N (Jonkers and Kučera, 2015) (Fig. 8a). With the reduction in (sea-)ice cover during the initial deglaciation, N. pachyderma starts occurring earlier in summer, persisting at the same low temperatures. As the deglaciation progresses, the “cold” population (P1), with a similar unimodal distribution as in the glacial, is joined by a second “warm” population (P2) that starts appearing in late summer. The isotopic difference between P1 and P2 (0.9±0.4 ‰) corresponds to a temperature offset of about ∼4C, the same as observed today at 59 N (Jonkers et al., 2013).

The modern seasonal succession of P1 and P2 generates the same bimodality we observe in the δ18O of the mixed N. pachyderma populations during the deglaciation in our core record (Fig. 7). Such bimodality may well be an expression of two genetically different but morphologically identical “cryptic species” among N. pachyderma (Bauch et al., 2003; Darling et al., 2000; Kučera and Darling, 2002). Indeed, morphologies are indistinguishable among our encrusted specimens from 250 to 300 µm both in our cored sediment and in modern N. pachyderma from the time series sediment traps at 59 N during both seasonal maxima, regardless of the size fractions used (Jonkers et al., 2010, 2013). Unfortunately, the lack of organic matter, due to the process of low temperature ashing, concentrate and isolate the mineral shells from organic matter leaving a clean residue for isotope and chemical analysis, limits the ability to genetically analyse trap specimens. At our core site, increasing temperatures would have first caused the disappearance of the “cold” water population (P1) (∼9.5 ka) followed shortly after by the disappearance of “warm” population (P2) (Fig. 7) when Holocene temperatures at this latitude exceed N. pachyderma's upper tolerance limit of approximately 10 C (Darling et al., 2006).

Figure 7Average δ18O value of single-shell populations for specimens of N. pachyderma and G. bulloides across the deglaciation and for the glacial interval. (a–b) Mean and standard deviation of distinct populations of N. pachyderma plotted against core depth. (c–d) Mean and standard deviation of distinct populations of G. bulloides plotted against core depth. Calculated values for P1 and P2, as determined from mixture analysis (Hammer et al., 2001). Vertical bars represent the standard deviation for each population; depths where multiple symbols are present are where it is not possible to distinguish statistically either one or more populations; these thus represent a single population of the sample to the left. Horizontal dashed lines represent the averages for P1 and P2; grey line is the total population average as would be reconstructed from pooled shell analysis. Above 300 cm (<10 ka), the Holocene, N. pachyderma has disappeared from the site.


Figure 8Seasonal succession in the modern North Atlantic. Top panels indicate fluxes of N. pachyderma (grey) and G. bulloides (blue) from sediment traps in (a) the polar Greenland–Norwegian Sea (Wolfteich, 1994), (b) subpolar Irminger Sea (Jonkers et al., 2010, 2013; Jonkers and Kučera, 2015) and (c) temperate mid-North Atlantic (Wolfteich, 1994); same labels as in Fig. 1. Middle and bottom panels represent the temperature and salinity from Ocean ReAnalysis System 4 (ORAS4) (Balmaseda et al., 2013). For panels (a) and (c), cups have been rearranged to progress from January to December. The fluxes, temperature and salinity given represent an averaging over the trap deployment period.


4.2 Alternative mechanisms and scenarios

Single-specimen isotope analysis permits unravelling of mixed sedimentary assemblages into their constituent components. Here, we show that the warming trend within the average δ18O of pooled N. pachyderma is directly caused by the emergence of a “warm” population (P2) shifting the mean isotopic value toward a warmer signal, concealing the continued existence of the original “cold” (P1) population. Within the northern North Atlantic, an abrupt change occurs from a single peak in production during the Last Glacial Maximum (LGM) to two populations that remain approximately 4 C apart throughout the deglaciation, inferring that the difference in δ18O is temperature driven, consistent with present-day observations from subpolar sediment trap time series. However, alternative scenarios that give the same or a similar solution for the existence of two populations can be envisaged. Below, we discuss other causal mechanisms that might be inferred from the data, including a low-salinity meltwater effect (Duplessy et al., 1991), bioturbation (Lougheed et al., 2018) and/or population dynamics (Mix, 1987; Roche et al., 2018).

4.2.1 Warming trend or meltwater pulse?

Reconstructions of the Δδ18Osw anomaly between the LGM and modern (Duplessy et al., 1991) suggest a series of regions above the southerly displaced Polar Front where freshwater and meltwater entered the North Atlantic in sufficient volumes to perturb the system, from continental ice meltwater and/or riverine input. Throughout the deglacial period, advances in the subtropical water masses and retreats of the Polar Front occurred. Repeated invasion of high-temperature and salinity waters into the Nordic Seas has shown that the deglacial period was inherently highly dynamic and thus unstable compared to the LGM as evidenced by isotopic (Duplessy et al., 1992; Kroon et al., 1997) and radiocarbon (Waelbroeck et al., 2001) measurements. Meltwater released into the northern North Atlantic during this time would have led to an increase in stratification and thus a decrease in sea surface temperature (SST) altering the evaporation–precipitation (E-P) balance that drives the poleward advection of subtropical water high in both temperature and salinity (Duplessy et al., 1992). The two populations found in our core during the deglaciation might have resulted from one seasonal population experiencing meltwater and a second seasonal population occurring before or after a meltwater event. The presence of continental IRD downcore in T88-3P, without a clear concomitant “spike” in the δ18O, referred to in the literature as a “meltwater spike” (Berger et al., 1977; Jones and Ruddiman, 1982) of either G. bulloides (Figs. 2 and 6) or N. pachyderma (Fig. 5), would suggest that the difference in δ18O between the two populations is dominated by temperature, consistent with previous studies showing no meltwater spike (Duplessy et al., 1996; Straub et al., 2013). Indeed, the presence of both foraminifera and IRD together downcore does not necessarily imply cohabitation of the same environment, as the modern seasonal maximum in polar shell productivity occurs prior to the arrival of meltwater from icebergs (Fig. 8). The extremely low values of continental ice (δ18O: −30 to −40 ‰) should lead to δ18O and salinity anomalies in surface waters, but sea ice formed from ocean water will have little impact on δ18O despite an impact on salinity. Therefore, a concordial meltwater δ18O signal and the presence of IRD is not compulsory (Duplessy et al., 1996) with increased sea-ice formation predicted to occur during periods of increased freshwater and extended Arctic Ocean area (Duplessy et al., 1996).

Figure 9Output of estimate of pooled specimen variance for T88-3P. Using the unmixed populations of N. pachyderma from the deglacial interval, based on single-shell measurements presented here (Fig. 7), a pooled synthetic measurement was created using the probability of each population and a synthesised normal distribution with the same mean and standard deviation. Estimates were made for (a) 5, (b) 10, (c) 20 and (d) 50 specimens per pooled analysis. For each sample, 10 000 replicates were produced, plotted here as a “heat map”; the colours represent the probability (counts normalised to 1) of a particular value occurring as a pooled value.


4.2.2 Spatial rather than temporal populations: shallow or deep?

Differences in depth habitat rather than timing might account for our observations. Depending on the structure of the water column, i.e. the depth of the surface mixed layer and the degree of stratification (see Metcalfe et al., 2015 for a discussion), the populations could represent one shallower and one deeper population that are not divided temporally but vertically within the water column (Fig. 8). Observations from the subpolar IRM time series sediment traps show that the first maximum occurs at an earlier time when the water column is well mixed, so that two vertically divided populations, i.e. one shallow and one deep, would have a similar δ18O signature. The second maximum in IRM occurs at a later time, during increased water column stratification; therefore, a shallow and deep population's δ18O should diverge. Therefore, only when the water column is stratified would it be possible to produce two theoretical populations different in δ18O, in much the same way as discrete species calcifying at different depths acquire an isotopic offset, enabling the use of Δδ18O as a proxy for past ocean stratification (Emiliani, 1954; Lototskaya and Ganssen, 1999; Mulitza et al., 1997). Following this line of reasoning, our results would suggest that the water column was more stratified during the deglaciation and well mixed during the LGM and Holocene.

One approach to further differentiate between depths is the carbon isotope (δ13C) signal, as seawater δ13C has a distinct signature, due in part to photosynthetic fractionation in the surface ocean enriching the euphotic zone in 13C. Exported organic matter may become remineralised at the base of the deep chlorophyll maximum enriching the euphotic zone in 12C at greater depth. Thus, the δ13C of foraminifera that have grown at different depths in these water masses should also have different values for each subpopulation, notwithstanding species specific vital effects. However, differences or similarities in carbon isotopes can also arise by alteration in food (either through grazing on different trophic levels and/or the types of food); if two seasonal populations of foraminifera existed, either their food source thrived for longer or a succession of more oligotrophic tolerant phytoplankton occurred. Therefore, the carbon isotope signature can be related to both scenarios.

By contrast, our results show no differences in δ13C signature between the two populations of N. pachyderma, notwithstanding interspecimen variance. What is directly observable, however, is that the IRM shows that there are two populations occurring seasonally. Second, there are no morphological differences observed between the specimens of N. pachyderma that isotopically belong to different populations in our core record, nor between the early and late summer maxima in N. pachyderma with a similarly distinct isotope composition in the modern sediment trap time series. Most species undergo wall thickening with depth (Brummer et al., 1986, 1987; Hemleben et al., 1985; Reynolds et al., 2018; Steinhardt et al., 2015) whilst some, including N. pachyderma, add a thick calcite crust with a different δ18O signature overprinting previous layers of the shell (Kozdon et al., 2009). This crust, however, is an ontogenetic feature (Brummer et al., 1986, 1987; Steinhardt et al., 2015) that is present in both seasons at IRM (Jonkers et al., 2010, 2013).

4.2.3 Sedimentary processes: dissolution and bioturbation

Seafloor processes such as dissolution and bioturbation may alter sediment populations in both isotope composition (Bard et al., 1987; Lougheed et al., 2017; Wit et al., 2013) and faunal composition (Bard, 2001; Löwemark, 2007; Löwemark et al., 2008). Dissolution not only removes “time” from the sediment but also leads to specimens being found together that have once been separated by centimetres of sedimentary material, as younger shells are deposited next to freshly exposed older shells (Lougheed et al., 2018). Similarly, depending on the oxygen content of sediments and the type and abundance of bottom fauna, bioturbation by benthic organisms can alter the sequence of cause and causality (Lougheed et al., 2018). Particle grain size distributions may also change due to bioturbation (Bard, 2001) if two species have differences in their absolute size, such as those measured here; it may show distinct isotope differences given species-specific size distributions (Brummer et al., 1986; Peeters et al., 1999). Thus, the two populations found in N. pachyderma δ18O could reflect relict specimens displaced in core depth, and therefore in time, given there is a shift in the sedimentation rate of core T88-3P between 290 and 410 cm and not one between 515 and 565 cm (Fig. 3). However, such sorting effects can be excluded here since both N. pachyderma populations and G. bulloides come from the same > 250 µm size fraction. It is important to note that for several depths in core this second population may only represent a few specimens (n<3specimens=3; and n<4specimens=5). Similarly, Löwemark et al. (2007, 2008) have shown that it is possible to have apparent differences due to the original abundance of the bioturbated species (Bard et al., 1987; Löwemark and Grootes, 2004). The difference between the two populations in G. bulloides and N. pachyderma could reflect a change in dominance of the foraminiferal assemblage, with bioturbation becoming more obvious in N. pachyderma as the species abundance reduces (Figs. 4 and 5).

Bioturbation is more easily detected during a period of pronounced climatic change, i.e. when the two end-member samples have the largest difference, as the signal of bioturbation becomes more pronounced. This may explain the difference in populations between the results from the deglacial (290 to 410 cm; Fig. 7) and glacial (515 to 565 cm; Fig. 7) sections despite both sections having similar abundance shifts (Figs. 4 to 6). As the resultant single-specimen δ18O distribution is a product of species-specific temperature tolerances (Mix, 1987; Roche et al., 2018), the visibility of bioturbation is especially enhanced at periods of sharp climatic transition. If the climatic signal crosses through a species temperature tolerance, then two separate warm and cold populations should exist separated both in time and core depth; bioturbation will then mix these populations together. However, we exclude this particular scenario because sedimentary features (Fig. 4) indicate a lack of discernible mixing, i.e. the sharpness of the IRD percentage and the XRF data indicate that bioturbation is at a minimum.

4.3 Palaeoceanographic implications: probability of drawing from either P1 or P2

The implications for the climate of the past are two-fold. First, our results suggest that there is more than one population of the polar N. pachyderma during the deglaciation and that its continued presence throughout much of this time period puts doubt on two discrete modes. The presence of both a colder population and warmer population suggests that this period is characterised by heightened seasonality, given that the climate conditions prevalent at ∼56 N supported two populations of N. pachyderma. This heightened “seasonality” is also visible in the increased standard deviation (Fig. 5) for these samples. The second implication is that this causal mechanism (i.e. seasonally distinct populations occurring during a climate transition) may not be captured using a pooled sample approach, given two distinct reactions to the same climate transition. It is important to note that G. bulloides also on occasion has more than a single population at this core site; however, the species' cosmopolitan and optimistic nature make it less surprising that expansion of seasonal variables that intersects the species' tolerances will lead to an expansion of its ecological range (e.g. Metcalfe et al., 2019). The appearance of a second population in N. pachyderma is more surprising because ecological expansion for this species can only be unidirectional (i.e. into water masses with higher temperatures) given its dominance of the polar environment. Therefore, given the use of N. pachyderma as a polar water mass indicator, we chose to investigate how multiple populations would impact pooled analysis. The unmixing algorithm used in this paper gives the probability of each distinct population, using each population and their calculated mean and standard deviation, to generate a normal distribution for the populations determined via statistical unmixing (Hammer et al., 2001; Wit et al., 2013). Using these data, it is possible to model the theoretical effect of sample size on the resultant stable isotope measurements (Morard et al., 2016). For simplicity, we assume that each specimen contributes an equal weighting to the overall pooled stable isotope value; of course, in reality, each specimen will contribute an amount of CO2 equal to its weight. This assumption will result in some error associated with our prediction of pooled specimen δ18O values due to kinetic fractionation during conversion from CaCO3 to CO2 (and H2O). The theoretical specimens were picked from either population or for those samples in which only a single population exists (at either limit of our sampling), using the rand function of MatLab. The function rand is statistically uniform throughout the range 0 to 1 and therefore can be used to construct a random number generator to define which population each theoretical specimen would have belonged to, using the following equations:

(2)r=randNpool,1pδ18Opop. II,(3)R1=normrndδ18Oμpop. I,δ18Oσpop. I,Npool,1,(4)R2=normrndδ18Oμpop. II,δ18Oσpop. II,Npool,1,(5)S=(1-r)R1+rR2.

The number of pooled specimens (in-group analysis; Npool) was varied between iterations of the model, so that 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90 and 100 draws/specimens were used for each subsequent iteration. For each depth, 10 000 redraws were performed, because we use a large number for resampling (N=10 000), and due to the central limit theorem, the average δ18O between the different iterations (variable number of pooled specimens) remains near constant. Therefore, the 2.5th, 25th, 75th and 97.5th quartiles were used to visually compare the spread of the data between different numbers of pooled specimens (Fig. S2 in the Supplement) and the probability of δ18O values occurring calculated (Fig. 9). The results of the model indicate that, for a foraminiferal fossil population composed of more than one discrete subpopulation, caution should be applied when using a small number of specimens for pooled analysis to ascertain an average state of the climate. Whilst the spread in values is narrower for some intervals, several downcore samples have a spread of 1 to 1.5 ‰ (Fig. 9).

The fact that a sample may have a single population for one or both species and intermittently two or more populations (i.e. N. pachyderma) may further complicate species comparison (e.g. Δδ18O). The emergence of a second population within N. pachyderma during the last deglaciation at the species southerly boundary indicates that other species with multiple abundance or size maxima (Schmidt et al., 2004a, b) may have a similarly hidden seasonal complexity within the stable isotope composition of pooled specimens. If these populations do not represent ecophenotypes, but instead are analogous to cryptic speciation in which populations are indistinguishable morphologically (Kučera and Darling, 2002; Morard et al., 2016), then pooled isotope measurement of such a sample will accidentally “pick” from multiple populations. Therefore, the wide use of N. pachyderma isotopes as a measure of sea-level rise, rate of deglaciation or ice-volume change based on the δ18O of pooled specimens may be unduly skewed.

5 Conclusions

Our findings expose and resolve the seasonal complexity that exists hidden in the δ18O produced within pooled specimens whilst highlighting the usefulness of integrating downcore studies with modern time series observation in the interpretation of species ecology for palaeoceanographic research. Using sediment trap time series data as modern keys to past climate conditions, our results imply that conditions existing today within the subpolar Irminger Sea prevailed at significantly more southerly latitudes throughout the last deglaciation. The remarkable difference between the transition (deglaciation) and the two climatic modes (glacial and interglacial) suggests that the mid-North Atlantic has an intermediate “deglacial” stable mode, rather than gradually shifting from glacial to interglacial. Our observation of a distinct bimodality throughout the deglaciation has important implications for how δ18O records can be interpreted given present-day seasonality; however, the interpretation of N. pachyderma as two populations instead of one is consistent with previous studies. Therefore, the common use of this species as a measure of sea-level rise, rate of deglaciation or ice-volume change, or ocean warming and stratification based on the δ18O of pooled specimens may be unduly skewed.

Data availability

The data produced from APNAP II T88-3P and used in this research can be accessed online (, Brummer et al., 2020).

Sample availability

Access to APNAP II T88-3P material is granted via request to Gerald Ganssen (VUA).


The supplement related to this article is available online at:

Author contributions

GMG was chief scientist of APNAP II (RV Tyro) during core retrieval, initiated and supervised the study, together with GJAB. WF conducted the investigation. Jv'tH performed abundance counts. GJAB, WF, BM and GMG contributed to data analysis and interpretation. BM made the figures and performed statistics. GJAB and BM wrote the manuscript with contributions from all authors.

Competing interests

The authors declare that they have no conflict of interest.


The captain and crew of the RV Tyro are thanked for core retrieval during APNAP II (1988). Core scanning was supported by the NWO programme on advanced instrumentation. Rineke Gieles is thanked for assisting with XRF core scanning. Sample treatment and preparation was conducted in the Sediment Laboratory of the Vrije Universiteit Amsterdam (VUA). Hubert Vonhof and Suzanne J. A. Verdegaal-Warmerdam (stable isotope laboratory of the VUA) are thanked for assistance with stable isotope analysis. Brett Metcalfe is supported by a Laboratoire d'excellence (LabEx) of the Institut Pierre-Simon Laplace (Labex L-IPSL). The authors wish to thank the Integrated Climate Data Center (ICDC, DE) for their online live access servers that provided access to atlas and reanalysis data used within this study.

Financial support

This research has been supported by the Netherlands Organisation for Scientific Research (SCAN2 (grant no. 834.11.003) and Digging for Density (grant no. 822.01.0.19)) and the Darwin Center for Biogeosciences (Sensing Seasonality (grant no. 3040)). Brett Metcalfe is supported by Laboratoire d'excellence (LabEx) of the Institut Pierre-Simon Laplace (LabEx L-IPSL), funded by the French Agence Nationale de la Recherche (grant no. ANR-10-LABX-0018).

Review statement

This paper was edited by Marit-Solveig Seidenkrantz and reviewed by two anonymous referees.


Balmaseda, M. A., Mogensen, K., and Weaver, A. T.: Evaluation of the ECMWF ocean reanalysis system ORAS4, Q. J. Roy. Meteorol. Soc., 139, 1132–1161, 2013. 

Bard, E.: Paleoceanographic implications of the difference in deep-sea sediment mixing between large and fine particles, Paleoceanography, 16, 235–239, 2001. 

Bard, E., Arnold, M., Duprat, J., Moyes, J., and Duplessy, J.-C.: Reconstruction of the last deglaciation: deconvolved records of δ18O profiles, micropaleontological variations and accelerator mass spectrometric 14C dating, Clim. Dynam., 1, 101–112, 1987. 

Bauch, D., Darling, K., Simstich, J., Bauch, H. A., Erlenkeuser, H., and Kroon, D.: Palaeoceanographic implications of genetic variation in living North Atlantic N. pachyderma, Nature, 424, 299–302, 2003. 

Berger, W. H.: Deep-sea carbonate and the deglaciation preservation spike in pteropods and foraminifera, Nature, 269, 301–304, 1977. 

Billups, K. and Spero, H. J.: Reconstructing the stable isotope geochemistry and paleotemperatures of the equatorial Atlantic during the last 150 000 years: Results from individual foraminifera, Paleoceanography, 11, 217–238,, 1996. 

Breitenbach, S. F. M. and Bernasconi, S. M.: Carbon and oxygen isotope analysis of small carbonate samples (20 to 100 µg) with a GasBench II preparation device, Rap. Commun. Mass Spec., 25, 1910–1914, 2011. 

Brummer, G.-J. A., Hemleben, C., and Spindler, M.: Planktonic foraminiferal ontogeny and new perspectives for micropalaeontology, Nature, 319, 50–52, 1986. 

Brummer, G.-J. A., Hemleben, C., and Spindler, M.: Ontogeny of extant spinose planktonic foraminifera (Globigerinidae): A concept exemplified by Globigerinoides sacculifer (Brady) and G. ruber (d'Orbigny), Mar. Micropaleontol., 12, 357–381, 1987. 

Brummer, G.-J. A., Metcalfe, B., Feldmeijer, W., Prins, M. A., van 't Hoff, J., and Ganssen, G. M.: Modal shift in North Atlantic seasonality during the last deglaciation (Version 1) [Data set], Climate of the Past, Zenodo,, 2020. 

Darling, K., Wade, C. M., Steward, I. A., Kroon, D., Dingle, R., and Leigh Brown, A. J.: Molecular evidence for genetic mixing of Arctic and Antarctic subpolar populations of planktonic foraminifers, Nature, 405, 43–47, 2000. 

Darling, K. F., Kučera, M., Kroon, D., and Wade, C. M.: A resolution for the coiling direction paradox in Neogloboquadrina pachyderma, Paleoceanography, 21, PA2011, 2006. 

Dempster, A. P., Laird, N. M., and Rubin, D. B.: Maximum likelihood from incomplete data via the EM algorithm, J. Roy. Stat. Soc. B, 39, 1–38, 1977. 

Duplessy, J. C., Labeyrie, L., Juilletleclerc, A., Maitre, F., Duprat, J., and Sarnthein, M.: Surface salinity reconstruction of the north-atlantic ocean during the last glacial maximum, Oceanologica Acta, 14, 311–324, 1991. 

Duplessy, J. C., Labeyrie, L., Arnold, M., Paterne, M., Duprat, J., and van Weering, T. C. E.: Changes in surface salinity of the North Atlantic Ocean during the last deglaciation, Nature, 358, 485–488, 1992. 

Duplessy, J. C., Labeyrie, L., and Paterne, M.: North Atlantic sea surface conditions during the Younger Dryas cold event, in: Late Quarternary Palaeoceanography of the North Atlantic Margins, edited by: Andrews, J. T., Austin, W. E. N., Bergsten, H., and Jennings, A. E., Geological Society Special Publication, 1996. 

Emiliani, C.: Depth habitats of some species of pelagic foraminifera as indicated by oxygen isotope ratios, Am. J. Sci., 252, 149–158, 1954. 

Epstein, S. and Mayeda, T.: Variation of O18 content of waters from natural sources, Geochim. Cosmochim. Acta, 4, 213–224, 1953. 

Falkowski, P. G. and Oliver, M. J.: Mix and match: how climate selects phytoplankton, Nat. Rev. Microbiol., 5, 813–819,, 2007. 

Feldmeijer, W., Metcalfe, B., Brummer, G. J. A., and Ganssen, G. M.: Reconstructing the depth of the permanent thermocline through the morphology and geochemistry of the deep dwelling planktonic foraminifer Globorotalia truncatulinoides, Paleoceanography, 30, 1–22, 2015. 

Frew, R. D., Dennis, P. F., Heywood, K. J., Meredith, M. P., and Boswell, S. M.: The oxygen isotope composition of water masses in the northern North Atlantic, Deep-Sea Res. I, 47, 2265–2286, 2000. 

Fritz-Endres, T., Dekens, P. S., Fehnrenbacher, J., Spero, H. J., and Stine, A.: Application of individual foraminifera Mg∕Ca and δ18O analyses for paleoceanographic reconstruction in active depositional environments, Paleoceanogr. Paleoclimatol., 34, 1610–1624,, 2019. 

Ganssen, G. M. and Kroon, D.: The isotopic signature of planktonic foraminifera from NE Atlantic surface sediments: implications for the reconstruction of past oceanic conditions, J. Geol. Soc., 157, 693–699, 2000. 

Ganssen, G. M., Peeters, F. J. C., Metcalfe, B., Anand, P., Jung, S. J. A., Kroon, D., and Brummer, G.-J. A.: Quantifying sea surface temperature ranges of the Arabian Sea for the past 20 000 years, Clim. Past, 7, 1337–1349,, 2011. 

Grant, K., Rohling, E., Bar-Matthews, M., Ayalon, A., Medina-Elizalde, M., Bronk Ramsey, C., Satow, C., and Roberts, A. P.: Rapid coupling between ice volume and polar temperature over the past 150,000 years, Nature, 491, 744–747,, 2012. 

Groeneveld, J., Ho, S. L., Mackensen, A., Mohtadi, M., and Laepple, T.: Deciphering the variability in Mg∕Ca and stable oxygen isotopesof individual foraminifera, Paleoceanogr. Paleoclimatol.,, 2019. 

Gruber, N., Keeling, C. D., and Bates, N. R.: Interannual Variability in the North Atlantic Ocean Carbon Sink, Science, 298, 2374–2378, 2002. 

Hammer, Ø., Harper, D. A. T., and Ryan, P. D.: PAST: Paleontological Statistics software package for education and data analysis, Palaeontologica Electronica, 4, 9 pp., 2001. 

Hemleben, C., Spindler, M., Breitinger, I., and Deuser, W. G.: Field and laboratory studies on the ontogeny and ecology of some globorotaliid species from the Sargasso Sea off Bermuda, J. Foramin. Res., 15, 254–272, 1985. 

Huybers, P.: Early Pleistocene Glacial Cycles and the Integrated Summer Insolation Forcing, Science, 313, 508–511, 2006. 

Jones, G. A. and Ruddiman, W. F.: Assessing the Global Meltwater Spike, Quarternary Res., 17, 148–172,, 1982. 

Jonkers, L. and Kučera, M.: Global analysis of seasonality in the shell flux of extant planktonic Foraminifera, Biogeosciences, 12, 2207–2226,, 2015. 

Jonkers, L., Brummer, G. J. A., Peeters, F. J. C., van Aken, H. M., and De Jong, M. F.: Seasonal stratification, shell flux, and oxygen isotope dynamics of left-coiling N. pachyderma and T. quinqueloba in the western subpolar North Atlantic, Paleoceanography, 25, PA2204, 2010. 

Jonkers, L., van Heuven, S., Zahn, R., and Peeters, F. J. C.: Seasonal patterns of shell flux, δ18O and δ13C of small and large N. pachyderma (s) and G. bulloides in the subpolar North Atlantic, Paleoceanography, 28, 164–174,, 2013. 

Killingley, J. S., Johnson, R. F., and Berger, W. H.: Oxygen and carbon isotopes of individual shells of planktonic foraminifera from Ontong-Java Plateau, Equatorial Pacific, Palaeogeogr. Palaeocl., 33, 193–204, 1981. 

Kohfeld, K. E., Fairbanks, R. G., Smith, S. L., and Walsh, I. D.: Neogloboquadrina pachyderma (sinistral coiling) as paleoceanographic tracers in polar oceans: Evidence from northeast water polynya plankton tows, sediment traps, and surface sediments, Paleoceanography, 11, 679–699, 1996,, 1996. 

Koutavas, A. and Joanides, S.: El Niño–Southern Oscillation extrema in the Holocene and Last Glacial Maximum, Paleoceanography, 27, PA4208,, 2012. 

Koutavas, A., deMenocal, P. B., Olive, G. C., and Lynch-Stieglitz, J.: Mid-Holocene El Niño–Southern Oscillation (ENSO) attenuation revealed by individual foraminifera in eastern tropical Pacific sediments, Geology, 34, 993,, 2006. 

Kozdon, R., Ushikubo, T., Kita, N. T., Spicuzza, M., and Valley, J. W.: Intratest oxygen isotope variability in the planktonic foraminifer N. pachyderma: Real vs. apparent vital effects by ion microprobe, Chem. Geol., 258, 327–337, 2009. 

Kretschmer, K., Kucera, M., and Schulz, M.: Modeling the distribution and seasonality of Neogloboquadrina pachyderma in the North Atlantic Ocean during Heinrich Stadial 1, Paleoceanography, 31, 986–1010, 2016. 

Kroon, D., Austin, W. E. N., Chapman, M. R., and Ganssen, G. M.: Deglacial surface circulation changes in the northeastern Atlantic: Temperature and salinity records off NW Scotland on a century scale, Paleoceanography, 12, 755–763, 1997. 

Kučera, M. and Darling, K. F.: Cryptic species of planktonic foraminifera: their effect on palaeoceanographic reconstructions, Philos. T. Roy. Soc. A, 360, 695–718, 2002. 

Kučera, M., Rosell-Melé, A., Schneider, R., Waelbroeck, C., and Weinelt, M.: Multiproxy approach for the reconstruction of the glacial ocean surface (MARGO), Quaternary Sci. Rev., 24, 813–819, 2005. 

Leduc, G., Vidal, L., Cartapanis, O., and Bard, E.: Modes of eastern equatorial Pacific thermocline variability: Implications for ENSO dynamics over the last glacial period, Paleoceanography, 24, PA3202,, 2009. 

Lototskaya, A. and Ganssen, G. M.: The structure of Termination II (penultimate deglaciation and Eemian) in the North Atlantic, Quaternary Sci. Rev., 18, 1641–1654, 1999. 

Lougheed, B. C. and Obrochta, S.: MatCal: Open Source Bayesian 14C Age Calibration in Matlab, J. Open Res. Softw., 4, e42,, 2016. 

Lougheed, B. C., Metcalfe, B., Ninnemann, U. S., and Wacker, L.: Moving beyond the age–depth model paradigm in deep-sea palaeoclimate archives: dual radiocarbon and stable isotope analysis on single foraminifera, Clim. Past, 14, 515–526,, 2018. 

Löwemark, L.: Importance and Usefulness of Trace fossils and Bioturbation in Paleoceanography, in: Trace Fossils: Concepts, Problems, Prospects, edited by: Miller, W., Elsevier, Amsterdam, 2007. 

Löwemark, L. and Grootes, P. M.: Large age differences between planktic foraminifers caused by abundance variations and Zoophycos bioturbation, Paleoceanography, 19, PA2001,, 2004. 

Löwemark, L., Konstantinou, K. I., and Steinke, S.: Bias in foraminiferal multispecies reconstructions of paleohydrographic conditions caused by foraminiferal abundance variations and bioturbational mixing: A model approach, Mar. Geol., 256, 101–106, 2008. 

Mikis, A., Hendry, K. R., Pike, J., Schmidt, D. N., Edgar, K. M., Peck, V., Peeters, F. J. C., Leng, M. J., Meredith, M. P., Todd, C. L., Stammerjohn, S., and Ducklow, H.: Temporal variability in foraminiferal morphology and geochemistry at the West Antarctic Peninsula: a sediment trap study, Biogeosciences, 16, 3267–3282,, 2019. 

Metcalfe, B., Feldmeijer, W., de Vringer-Picon, M., Brummer, G.-J. A., Peeters, F. J. C., and Ganssen, G. M.: Late Pleistocene glacial–interglacial shell-size–isotope variability in planktonic foraminifera as a function of local hydrography, Biogeosciences, 12, 4781–4807,, 2015. 

Metcalfe, B., Feldmeijer, W., and Ganssen, G. M.: Oxygen isotope variability of planktonic foraminifera provide clues to past upper ocean seasonal variability, Paleoceanography and Paleoclimatology, 34, 374–393,, 2019. 

Mix, A. C.: The oxygen-isotope record of deglaciation, in: North America and adjacent oceans during the last deglaciation, in: The Geology of America, edited by: Ruddiman, W. F. and Wright, H. E. J., Geological Society of America, Boulder, Colorado, 1987. 

Morard, R., Reinelt, M., Chiessi, C. M., Groeneveld, J., and Kucera, M.: Tracing shifts of oceanic fronts using the cryptic diversity of the planktonic foraminifera Globorotalia inflata, Paleoceanography, 31, 1193–1205, 2016. 

Mulitza, S., Dürkoop, A., Hale, W., Wefer, G., and Niebler, H. S.: Planktonic foraminifera as recorders of past surface-water stratification, Geology, 25, 335–338, 1997. 

Oba, T.: Paleoceanographic information obtained by the isotopic measurement of individual foraminiferal specimens, in: Proceedings First International Conference Asian Marine Geology, Shanghai, 1988, China Ocean Press, Beijing, 169–180, 1990. 

Oba, T.: Oxygen and carbon isotopic composition of planktonic foraminifera tests collected with sediment traps from the Japan Trench, La mer – Societe franco-japonaise d’oceanographie, 29, 190–192, 1991. 

Pearson, P.: Oxygen isotopes in foraminifera: Overview and historical review, Reconstructing Earth's Deep-Time Climate – The State of the Art in 2012, Paleontological Society Short Course, 1–38, 2012. 

Peeters, F. J. C., Ivanova, E., Conan, S. M. H., Brummer, G.-J. A., Ganssen, G. M., Troelstra, S., and van Hinte, J.: A size analysis of planktic foraminifera from the Arabian Sea, Mar. Micropaleontol., 36, 31–63, 1999. 

Pracht, H., Metcalfe, B., and Peeters, F. J. C.: Oxygen isotope composition of the final chamber of planktic foraminifera provides evidence of vertical migration and depth-integrated growth, Biogeosciences, 16, 643–661,, 2019. 

Rasmussen, S. O., Seierstad, I. K., Andersen, K. K., Bigler, M., Dahl-Jensen, D., and Johnsen, S. J.: Synchronization of the NGRIP, GRIP, and GISP2 ice cores across MIS 2 and palaeoclimatic implications, Quaternary Sci. Rev., 27, 18–28, 2008. 

Reimer, P. J., Bard, E., Bayliss, A., Beck, J. W., Blackwell, P. G., Bronk Ramsey, C., Buck, C. E., Cheng, H., Edwards, R. L., and Friedrich, M.: IntCal13 and Marine13 Radiocarbon age calibration curves 0–50,000 years cal BP, Radiocarbon, 55, 1869–1887, 2013. 

Reynolds, C. E., Richey, J. N., Fehrenbacher, J. S., Rosenheim, B. E., and Spero, H. J.: Environmental controls on the geochemistry of Globorotalia truncatulinoides in the Gulf of Mexico: Implications for paleoceanographic reconstructions, Mar. Micropaleontol., 142, 92–104, 2018. 

Richter, T. O., van der Gaast, S., Koster, B., Vaars, A., Gieles, R., de Stigter, H. C., De Haas, H., and van Weering, T. C. E.: The Avaatech XRF Core Scanner: technical description and applications to NE Atlantic sediments, Geol. Soc. London Spec. Publ., 267, 39–50, 2006. 

Roche, D. M., Waelbroeck, C., Metcalfe, B., and Caley, T.: FAME (v1.0): a simple module to simulate the effect of planktonic foraminifer species-specific habitat on their oxygen isotopic content, Geosci. Model Dev., 11, 3587–3603,, 2018. 

Schiffelbein, P. and Hills, S.: Direct assessment of stable isotope variability in planktonic foraminifera populations, Palaeogeography, Palaeoclimatology, Palaeoecology, 48, 197–213,, 1984. 

Schmidt, D. N., Renaud, S., Bollmann, J., Schiebel, R., and Thierstein, H. R.: Size distribution of Holocene planktic foraminifer assemblages: biogeography, ecology and adaptation, Mar. Micropaleontol., 50, 319–338, 2004a. 

Schmidt, D. N., Thierstein, H. R., Bollmann, J., and Schiebel, R.: Abiotic Forcing of Plankton Evolution in the Cenozoic, Science, 303, 207–210, 2004b. 

Schmittner, A., Urban, N. M., Shakun, J. D., Mahowald, N. M., Clark, P. U., Bartlein, P. J., Mix, A. C., and Rosell-Melé, A.: Climate sensitivity estimated from temperature reconstructions of the Last Glacial Maximum, Science, 334, 1385–1388, 2011. 

Scussolini, P., van Sebille, E., and Durgadoo, J. V.: Paleo Agulhas rings enter the subtropical gyre during the penultimate deglaciation, Clim. Past, 9, 2631–2639,, 2013. 

Shackleton, N. J.: Oxygen isotopes, ice volume and sea level, Quaternary Sci. Rev., 6, 183–190,, 1987. 

Spötl, C. and Vennemann, T. W.: Continuous-flow isotope ratio mass spectrometric analysis of carbonate minerals, Rapid Commun. Mass Spectrom., 17, 1004–1006, 2003. 

Steinhardt, J., de Nooijer, L., Brummer, G.-J. A., and Reichart, G. J.: Profiling planktonic foraminiferal crust formation, Geochem. Geophys. Geosy., 16, 2409–2430, 2015. 

Straub, M., Tremblay, M. M., Sigman, D. M., Studer, A. S., Ren, H., Toggweiler, J. R., and Haug, G. H.: Nutrient conditions in the subpolar North Atlantic during the last glacial period reconstructed from foraminifera-bound nitrogen isotopes, Paleoceanography, 28, 79–90, 2013. 

Takagi, H., Moriya, K., Ishimura, T., Suzuki, A., Kawahata, H., and Hirano, H.: Exploring photosymbiotic ecology of planktc foraminifers from chamber-by-chamber isotopic history of individual foraminifers, Paleobiology, 41, 108–121, 2015. 

Takagi, H., Moriya, K., Ishimura, T., Suzuki, A., Kawahata, H., and Hirano, H.: Individual migration pathways of modern planktic foraminifers: Chamber-by-chamber assessment of stable isotopes, Paleontol. Res., 20, 268–284, 2016. 

van Sebille, E., Scussolini, P., Durgadoo, J. V., Peeters, F. J. C., Biastoch, A., Weijer, W., Turney, C., Paris, C. B., and Zahn R.: Ocean currents generate large footprints in marine palaeoclimate proxies, Nat. Commun., 6, 6521,, 2015.  

Waelbroeck, C., Duplessy, J.-C., Michel, E., Labeyrie, L., Paillard, D., and Duprat, J.: The timing of the last deglaciation in North Atlantic climate records, Nature, 412, 724–727, 2001. 

Waelbroeck, C., Labeyrie, L., Michel, E. , Duplessy, J. C., McManus, J. F., Lambeck, K., Balbon, E., and Labracherie, M.: Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records, Quaternary Sci. Rev., 21, 295–305,, 2002. 

Waelbroeck, C., Mulitza, S., Spero, H., Dokken, T., Kiefer, T., and Cortijo, E.: A global compilation of late Holocene planktonic foraminiferal 18O: relationship between surface water temperature and 18O, Quaternary Sci. Rev., 24, 853–868, 2005. 

Weltje, G. J.: End-member modeling of compositional data: Numerical-statistical algorithms for solving the explicit mixing problem, Mathemat. Geol., 29, 503–549, 1997. 

Weltje, G. J. and Prins, M. A.: Muddled or mixed? Inferring palaeoclimate from size distributions of deep-sea clastics, Sediment. Geol., 162, 39–62, 2003. 

Weltje, G. J. and Tjallingii, R.: Calibration of XRF core scanners for quantitative geochemical logging of sediment cores: Theory and application, Earth Planet. Sc. Lett., 274, 423–438, 2008. 

Wit, J. C., Reichart, G. J. A., Jung, S. J., and Kroon, D.: Approaches to unravel seasonality in sea surface temperatures using paired single-specimen foraminiferal δ18O and Mg∕Ca analyses, Paleoceanography, 25, PA4220,, 2010. 

Wit, J. C., Reichart, G. J., and Ganssen, G. M.: Unmixing of stable isotope signals using single specimen δ18O analyses, Geochem. Geophys. Geosy., 14, 1312–1320, 2013. 

Wolfteich, C. M.: Satellite-derived sea surface temperature, mesoscale variability, and foraminiferal production in the North Atlantic, M.Sc., Massachusetts Institute of Technology/Woodshole Oceanographic Institution, 1994, 80 pp., 1994. 

Short summary
Here, mid-ocean seasonality is resolved through time, using differences in the oxygen isotope composition between individual shells of the commonly used (sub)polar planktonic foraminifera species in ocean-climate reconstruction, N. pachyderma and G. bulloides. Single-specimen isotope measurements during the deglacial period revealed a surprising bimodality, the cause of which was investigated.