Modal shift in North Atlantic seasonality during the last deglaciation

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 (δ¹⁸O) 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 δ¹⁸O of polar N. pachyderma reveals a steeper rate of ocean warming during the last deglaciation than appears from conventional pooled δ¹⁸O average values.

1 Introduction

1.1 Seasonality and single-foraminiferal analysis (SFA)

Stable oxygen isotopes (δ¹⁸O) 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 1 Map 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 ∼75 ^∘C 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

$$\begin{array}{}\text{(1)}& \text{Ratio of NPS}={\displaystyle \frac{\mathit{\text{N. pachyderma}}}{(\mathit{\text{N. pachyderma}}+\mathit{\text{G. bulloides}})}}.\end{array}$$

2.3 Single-foraminifera stable isotope geochemistry (δ¹⁸O; δ¹³C)

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 H₃PO₄ (45 ^∘C for 160 min). The resultant CO₂–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 CO₂ 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 (δ¹⁸O) 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 δ¹⁸O and δ¹³C, 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: n_total=414; n_{outlier corrected}=388; n_outlier=26; G. bulloides: n_total=439; n_{outlier corrected}=424; n_outlier=15) and glacial sections (N. pachyderma: n_total=490; n_{outlier corrected}=474; n_outlier=16; G. bulloides: n_total=496; n_{outlier corrected}=474; n_outlier=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 2 Raw 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 3 Age 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 δ¹⁸O tie points.

Figure 4 Sediment 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 5 Single-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 6 Single-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 ¹⁴C yr, with an error of 200 ¹⁴C yr, expressed mathematically as ΔR: 0±200 ¹⁴C 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 1 Raw and calibrated radiocarbon ages. Conventional radiocarbon age represents the measured radiocarbon age (see Table S1 in the Supplement) corrected for isotopic fraction.

2.4.2 Age model construction

Two age models were produced that used two independent methods: first, a radiocarbon-only age model, and second, δ¹⁸O 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⁻¹) to a rapid interglacial SAR (∼30 cm kyr⁻¹); 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 δ¹⁸O 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 δ¹⁸O and δ¹³C 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 δ¹⁸O 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 δ¹⁸O 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⁻¹, 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 δ¹⁸O 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 δ¹⁸O 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 δ¹⁸O: 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 δ¹⁸O average values and standard deviation lie between 2.5 and 4.5 ‰, with a remarkable consistent spread in the values; the lightest δ¹⁸O 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 δ¹⁸O average values and standard deviation during this interval range from 5.5 to 2.5 ‰.

3.3 Single-shell δ¹⁸O: 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 δ¹⁸O 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 δ¹⁸O average values and standard deviation during the deglacial interval range from 4 to 1.5 ‰.

3.4 Single-shell δ¹⁸O 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; n_groups=26; n_{within group}=474) and uncorrected for outliers (μ=0.55; min = 0.28; max = 0.85; σ=0.15; n_groups=26; n_{within group}=490) lower than the deglacial samples (300–420 cm) either uncorrected for outliers (μ=0.70; min = 0.26; max = 1.52; σ=0.31; n_groups=22; n_{within group}=414) or corrected for outliers (μ=0.49; min = 0.13; max = 1.11; σ=0.26; n_groups=22; n_{within 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; n_groups=23; n_{within group}=439) and corrected (μ=0.38; min = 0.19; max = 0.65; σ=0.13; n_groups=23; n_{within group}=424) are somewhat similar to the glacial uncorrected (μ=0.69; min = 0.33; max = 1.53; σ=0.28; n_groups=26; n_{within group}=496) and corrected (μ=0.48; min = 0.24; max = 0.88; σ=0.16; n_groups=26; n_{within group}=474) data.

3.5 Single-shell δ¹⁸O populations

Our results show that δ¹⁸O 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 δ¹⁸O values remains predominately unimodal throughout the deglacial section (13 out of 23 samples), the δ¹⁸O 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 δ¹⁸O persisting from the glacial (population P1) and a second population low in δ¹⁸O appearing at the onset of the deglaciation (population P2). The difference in δ¹⁸O 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 δ¹⁸O values of N. pachyderma become once more unimodal, now for P2, shortly before disappearing entirely until the present day. Carbon isotope values (δ¹³C) measured on the same shells of N. pachyderma do not appear to show this bimodal distribution, precluding the possibility that the two populations in δ¹⁸O represent a similar season but grew their shells at different depths given the enrichment and depletion with depth in seawater ¹³C associated with phytoplankton growth and decay. Since the δ¹⁸O 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 ∼4 ^∘C, 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 δ¹⁸O 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 7 Average δ¹⁸O 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 ($<\sim \mathrm{10}$ ka), the Holocene, N. pachyderma has disappeared from the site.

Figure 8 Seasonal 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 δ¹⁸O 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 δ¹⁸O 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 $\mathrm{\Delta }{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{sw}}$ 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 δ¹⁸O, 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 δ¹⁸O 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 (δ¹⁸O: −30 to −40 ‰) should lead to δ¹⁸O and salinity anomalies in surface waters, but sea ice formed from ocean water will have little impact on δ¹⁸O despite an impact on salinity. Therefore, a concordial meltwater δ¹⁸O 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 9 Output 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 δ¹⁸O signature. The second maximum in IRM occurs at a later time, during increased water column stratification; therefore, a shallow and deep population's δ¹⁸O should diverge. Therefore, only when the water column is stratified would it be possible to produce two theoretical populations different in δ¹⁸O, in much the same way as discrete species calcifying at different depths acquire an isotopic offset, enabling the use of Δδ¹⁸O 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 (δ¹³C) signal, as seawater δ¹³C has a distinct signature, due in part to photosynthetic fractionation in the surface ocean enriching the euphotic zone in ¹³C. Exported organic matter may become remineralised at the base of the deep chlorophyll maximum enriching the euphotic zone in ¹²C at greater depth. Thus, the δ¹³C 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 δ¹³C 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 δ¹⁸O 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 δ¹⁸O 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_{<\mathrm{3}\text{specimens}}=\mathrm{3}$; and $n_{<\mathrm{4}\text{specimens}}=\mathrm{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 δ¹⁸O 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 CO₂ equal to its weight. This assumption will result in some error associated with our prediction of pooled specimen δ¹⁸O values due to kinetic fractionation during conversion from CaCO₃ to CO₂ (and H₂O). 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:

$$\begin{array}{}\text{(2)}& {\displaystyle }r=\mathrm{rand}\left(N_{\mathrm{pool}},\mathrm{1}\right)\ge p\left({{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\text{pop. II}}\right),\text{(3)}& {\displaystyle }R\mathrm{1}=\mathrm{normrnd}\left({{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathit{\mu }\text{pop. I}},{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathit{\sigma }\text{pop. I}},N_{\mathrm{pool}},\mathrm{1}\right),\text{(4)}& {\displaystyle }R\mathrm{2}=\mathrm{normrnd}\left({{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathit{\mu }\text{pop. II}},{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathit{\sigma }\text{pop. II}},N_{\mathrm{pool}},\mathrm{1}\right),\text{(5)}& {\displaystyle }S=(\mathrm{1}-r)\cdot R\mathrm{1}+r\cdot R\mathrm{2}.\end{array}$$

The number of pooled specimens (in-group analysis; N_pool) 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 δ¹⁸O 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 δ¹⁸O 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. Δδ¹⁸O). 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 δ¹⁸O of pooled specimens may be unduly skewed.

5 Conclusions

Our findings expose and resolve the seasonal complexity that exists hidden in the δ¹⁸O 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 δ¹⁸O 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 δ¹⁸O of pooled specimens may be unduly skewed.