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


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 leftcoiling N. pachyderma (e.g. Mikis 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.

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.

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

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 Ratio of NPS = N. pachyderma (N. pachyderma + G. bulloides) .
2.3 Single-foraminifera stable isotope geochemistry (δ 18 O; δ 13 C) For single-shell stable isotope analysis, a continuous-flow isotope ratio mass spectrometer was used 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 3 PO 4 (45 • C for 160 min). The resultant CO 2 -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 2 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 (δ 18 O) in per mil (‰) versus Vienna Pee Dee belemnite (V-PDB). The reproducibility of an international car-bonate standard (IAEA-CO1) analysed is < 0.12 ‰ (1σ ) for both δ 18 O and δ 13 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 outliercorrected 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.
2.4 Core stratigraphy and age model

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  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 Ma-rine13 calibration curve (Reimer et al., 2013) and a reservoir age of 400 14 C yr, with an error of 200 14 C yr, expressed mathematically as R: 0 ± 200 14 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.

Age model construction
Two age models were produced that used two independent methods: first, a radiocarbon-only age model, and second, δ 18 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 −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 us-ing δ 18 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 δ 18 O and δ 13 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 δ 18 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 δ 18 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 di-  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 δ 18 O tie points. 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. verge. 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.

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 δ 18 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 outliercorrected 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.

Modern sediment trap record and ocean reanalysis
As the modern analogue of our distinct isotopic "endmembers", 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., 2010Jonkers et al., , 2013Jonkers 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.

IRD, abundance
The upper ∼ 290 cm of core T88-3P is Holocene in age as evidenced by near-uniform values of pooled specimen δ 18 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.

Single-shell δ 18 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 δ 18 O average values and standard deviation lie between 2.5 and 4.5 ‰, with a remarkable consistent spread in the values; the lightest δ 18 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 δ 18 O average values and standard deviation during this interval range from 5.5 to 2.5 ‰.

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

Single-shell δ 18 O populations
Our results show that δ 18 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 δ 18 O values remains predominately unimodal throughout the deglacial section (13 out of 23 samples), the δ 18 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 δ 18 O persisting from the glacial (population P1) and a second population low in δ 18 O appearing at the onset of the deglaciation (population P2). The difference in δ 18 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 δ 18 O values of N. pachyderma become once more unimodal, now for P2, shortly before disappearing entirely until the present day. Carbon isotope values (δ 13 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 δ 18 O represent a similar season but grew their shells at different depths given the enrichment and depletion with depth in seawater 13 C associated with phytoplankton growth and decay. Since the δ 18 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).

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(Jonkers et al., , 2013 (Fig. 8b). Neither of the N. pachyderma populations from IRM displays significant morphological differences (Jonkers et al., 2010(Jonkers et al., , 2013. By contrast, modern shell fluxes in the temperate North Atlantic at 48 • N, close to our core site, are domi-nated 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 δ 18 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(Jonkers et al., , 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).

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 δ 18 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 δ 18 O is temperature driven, consistent with presentday 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).

Warming trend or meltwater pulse?
Reconstructions of the δ 18 O 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 evaporationprecipitation (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 δ 18 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 δ 18 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 (δ 18 O: −30 to −40 ‰) should lead to δ 18 O and salinity anomalies in surface waters, but sea ice formed from ocean water will have little impact on δ 18 O despite an impact on salinity. Therefore, a concordial meltwater δ 18 O signal and the presence of IRD is not compulsory (Duplessy et al., 1996) with increased sea-ice formation predicted to  (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. occur during periods of increased freshwater and extended Arctic Ocean area (Duplessy et al., 1996).

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 δ 18 O signature. The second maximum in IRM occurs at a later time, during increased water column stratification; therefore, a shal-low and deep population's δ 18 O should diverge. Therefore, only when the water column is stratified would it be possible to produce two theoretical populations different in δ 18 O, in much the same way as discrete species calcifying at different depths acquire an isotopic offset, enabling the use of δ 18 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 (δ 13 C) signal, as seawater δ 13 C has a distinct signature, due in part to photosynthetic fractionation in the surface ocean enriching the euphotic zone in 13 C. Exported organic matter may become remineralised at the base of the deep chlorophyll maximum enriching the euphotic zone in 12 C at greater depth. Thus, the δ 13 C of foraminifera that have grown at different depths in these water masses should also have different values for each subpopulation, notwithwww.clim-past.net/16/265/2020/ Clim. Past, 16, 265-282, 2020  (Jonkers et al., 2010(Jonkers et al., , 2013Jonkers 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.
standing 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 δ 13 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(Brummer et al., , 1987Hemleben et al., 1985;Reynolds et al., 2018;Steinhardt et al., 2015) whilst some, including N. pachyderma, add a thick calcite crust with a different δ 18 O signature overprinting previous layers of the shell (Kozdon et al., 2009). This crust, however, is an ontogenetic feature (Brummer et al., 1986(Brummer et al., , 1987Steinhardt et al., 2015) that is present in both seasons at IRM (Jonkers et al., 2010(Jonkers et al., , 2013.

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 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.
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 δ 18 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 <3 specimens = 3; and n <4 specimens = 5). Similarly, Löwemark et al. (2007Löwemark et al. ( , 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 δ 18 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. 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 2 equal to its weight. This assumption will result in some error associated with our prediction of pooled specimen δ 18 O values due to kinetic fractionation during conversion from CaCO 3 to CO 2 (and H 2 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: 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 δ 18 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 δ 18 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. δ 18 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 δ 18 O of pooled specimens may be unduly skewed.

Conclusions
Our findings expose and resolve the seasonal complexity that exists hidden in the δ 18 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 δ 18 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 δ 18 O of pooled specimens may be unduly skewed.
Sample availability. Access to APNAP II T88-3P material is granted via request to Gerald Ganssen (VUA).
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.