Morphotypical and Geochemical Variations of Planktic Foraminiferal Species in Siberian and Central Arctic Ocean Core Tops Open Access

Polar amplification dictates that Arctic air temperatures will warm faster relative to lower latitudes in response to anthropogenic greenhouse warming. Coupled ocean-atmospheric general circulation models typically show a 2–3x increase (Holland & Bitz, 2003), with new studies showing a 4x average increase (Rantanen et al., 2022), in polar warming relative to the global average at twice the preindustrial atmospheric CO₂ level, which jeopardizes Arctic sea and land ice. Sea ice loss is predicted to be spatially variable and affect atmospheric climate modes such as the North Atlantic Oscillation (Holland & Bitz, 2003; Vihma, 2014).

In addition, near-future Greenland Ice Sheet (GrIS) stability is debated, with concerns that the GrIS may lose mass more quickly than predicted in response to global warming (Schaefer et al., 2016; Rantanen et al., 2022). The impacts of freshwater forcing, from melting Arctic land ice on the Atlantic Meridional Overturning Circulation (AMOC) is actively investigated, but short-term changes in Arctic oceanography could impact ocean circulation across longer timescales; increased freshwater from land ice melt in the Arctic is projected to slow mixing in the North Atlantic (Briner et al., 2020; He & Clark, 2022; Lofverstrom et al., 2022). Studies have also predicted a loss of Arctic ice analogous to early Holocene glacial melt and across other periods of ice sheet instability and suggest that AMOC could be significantly reduced (He & Clark, 2022; Lofverstrom et al., 2022). Therefore, it is imperative we understand how Arctic oceanography, regional land and sea ice distributions, and ocean temperatures have changed in response to past climatic forcing.

Many studies have used planktic foraminifera found in marine sediments as indicators of past climate change in the Arctic Ocean. Planktic foraminifera are ubiquitous, microscopic, marine protists that inhabit upper oceanic waters and build calcareous tests (shells), whose geochemistry and relative species abundances have been the backbone of paleoceanographic reconstructions over the past seventy years. Planktic foraminifera show a large latitudinal diversity gradient (Ruddiman, 1969; Schiebel & Hemleben, 2005), and species diversity decreases at the poles. Neogloboquadrina pachyderma (Ehrenberg, 1862) is the dominant planktic species found in polar latitudes, and as such, most high-latitude records of foraminiferal geochemistry typically focus on this species. The species is tolerant of cold temperatures (as low as –3°C) and a range of salinities from 30–69 (Hilbrecht, 1996; Bertlich et al., 2021), and has also been found at subtropical latitudes (Kennett & Srinivasan, 1983). These qualities establish N. pachyderma as a widely used indicator of Arctic paleoclimatic change (Kohfeld et al., 1996; Hillaire-Marcel et al., 2004; Risebrobakken & Berben, 2018).

The freezing water temperatures in the Arctic, low productivity, and large salinity gradients create a suboptimal environment for the less common Arctic planktic species Neogloboquadrina incompta (Parker, 1962), Turborotalita quinqueloba (Natland, 1938), and Turborotalita humilis (Brady, 1884; Bond et al., 2001; Darling et al., 2006; Cronin et al., 2014). Therefore, the abundance of these taxa, alongside N. pachyderma, are used as a proxy for changes in water mass circulation in the high-latitude oceans (Simstich et al., 2003). Neogloboquadrina incompta and N. pachyderma appear superficially similar in morphology (except for their coiling direction) but are genetically distinct (Darling et al., 2006; Schiebel & Hemleben, 2017). It is important to note that both species can have a minor proportion of individuals within the population that coil in the opposite direction to their standard orientation (Darling et al., 2006). Whereas N. pachyderma is typically sinistral (left-coiling), it can have a limited occurrence of ∼ 3% dextral (right-coiling) individuals (Darling et al., 2006). Additionally, there is evidence of higher percentages of aberrant coiling direction in laboratory populations, however it is not currently understood what implications this has on paleontological or modern observations (Davis et al., 2020).

Neogloboquadrina pachyderma abundances and isotopic geochemistry have been shown to reflect upper ocean conditions at marginal Arctic Ocean sites. Plankton tows in the Greenland Sea and Nansen Basin indicate that N. pachyderma size and geochemistry change with depth, and that the stable oxygen isotope composition (δ¹⁸O) of N. pachyderma shells is controlled by the temperature and the seawater δ¹⁸O (hereafter δ¹⁸O_seawater) of the water mass in which they calcify (Bauch et al., 1997; Kohfeld et al., 1996). Neogloboquadrina pachyderma δ¹⁸O consistently shows an offset of ∼ 1‰ from predicted equilibrium calcite δ¹⁸O, which is attributed to ‘vital effects’ (Bauch et al., 1997; Kohfeld et al., 1996). However, a recent study on small groups (four per measurement) of sub-Arctic foraminifera suggests a relatively minor degree of δ¹⁸O variability (0.11‰) unexplained by temperature and δ¹⁸O_seawater (Jonkers et al., 2022). In addition to geochemistry, maximum abundances (relative to depth at a local site) of living N. pachyderma follow modern chlorophyll maxima in the water column and have been observed with green cytoplasm, which suggests they consume primary producers in the Arctic (Kohfeld et al., 1996; Bauch et al., 1997). Recent work has identified diatoms as their primary food source in the Arctic, with N. pachyderma consuming various species based on availability (Greco et al., 2021).

Relative abundances of N. pachyderma and Turborotalita quinqueloba have also been proposed to track changes in water mass mixing over geologic time (Risebrobakken & Berben, 2018). Furthermore, N. pachyderma δ¹⁸O excursions have been used to define large surface freshwater fluxes in the Arctic Ocean in the Holocene (Poore et al., 1999; Bond et al., 2001). Reconstructions using N. pachyderma typically assume stationarity in its calcification preferences. However, the geochemical and morphotypical variability occurring among various N. pachyderma specimens and their relationship to calcification depths has not yet been fully explored across the Arctic (Altuna et al., 2018).

Altuna et al. (2018) identified five morphotypes within N. pachyderma (termed Nps-1 through Nps-5, following their convention). Previous work has also delineated N. pachyderma forms solely into “encrusted” and “nonencrusted” categories (Kohfeld et al., 1996) based on the putative presence of a thick calcite crust. “Encrusted” forms have been repeatedly observed deeper in the water column than “nonencrusted” forms, pointing toward vertical migration as an influence on N. pachyderma calcification (Kohfeld et al., 1996; Hillaire-Marcel et al., 2004; Altuna et al., 2018; Tell et al., 2022). While the species itself is often attributed to a range of depth habitats from 0 to 200 m, with maximum abundances around 50–100 m, various morphotypes are hypothesized to live at different depths within that range (Bauch et al., 1997; Volkmann, 2000; Bond et al., 2001; Altuna et al., 2018). Furthermore, the question of depth habitat as it relates to morphotype size has two conflicting hypotheses: one suggests larger forms live deeper in the water column, around 100–200 m, (Hillaire-Marcel et al., 2004), and the other suggests the inverse (Xiao et al., 2014). If morphotypes of N. pachyderma indeed calcify at different depths of the water column, they would exhibit different δ¹⁸O values relative to changes in temperature and δ¹⁸O_seawater within the water column. Vertical changes are steep in the Arctic Ocean but exhibit spatial variability. Generally, there is a strong halocline between a low-salinity mixed layer influenced by sea ice melt and riverine input versus higher-salinity Atlantic waters. The depth of the halocline would influence N. pachyderma δ¹⁸O. Thus, it is important to resolve the depths at which N. pachyderma and its morphotypes live and calcify across the Arctic Ocean to accurately interpret downcore records. While there is evidence for a relationship between N. pachyderma morphotypes, size, and stable isotopic composition in the Canadian Archipelago (Altuna et al., 2018), that environment does not reflect the entirety of Arctic oceanography. Here, we investigate multiple planktic species across mid- to late Holocene sediment obtained from core tops from the eastern and central Arctic Ocean and evaluate the implications of their geochemistry for calcification depth habitat.

For this study, we used nine core top sediments (0–1 cm and 1–2 cm; Table 1) spanning 76–90°N and 144°W–40°E, along an approximate transect from the Bering Strait toward the Fram Strait via the central Arctic Ocean (Fig. 1, Fig. A1). We obtained temperature and salinity data from the World Ocean Atlas (WOA; Locarnini et al., 2019; Zweng et al., 2019) for the upper 200 m of the water column (Fig. 1, Figs. A2–A4). World Ocean Atlas data are sourced from CTD casts and observations are interpolated to 102 standard depths. The dataset spans 1955 to present with a monthly resolution and has a maximum square grid resolution of 0.25°. In addition to WOA data, we collated CTD data from seven sites close to our core sites (Anderson, 2006; Swift, 2006; Rabe & Wisotzki, 2010; Björk, 2017). We then imaged and analyzed individual tests (Fig. 2), prepared assemblages for species and morphotype abundance (Fig. 3), and calculated benthic-to-planktic abundance ratios. Subsequently, we analyzed the stable isotope geochemistry of multiple neogloboquadrinids (Figs. 4–6), and Mg/Ca in N. pachyderma as proxies for temperature and δ¹⁸O_seawater (Fig. 7). Bulk dry samples were obtained from the United States Geologic Survey, Florence Bascom Geosciences Center (Table 1). Cores were originally collected during the following cruises: Polarstern Ark-VIII/3, the 1994 Arctic Ocean Section, and the Swedish-Russian-US Arctic Ocean Investigation of Climate-Cryosphere-Carbon Interactions (SWERUS) Leg 2 (The participants of ARK-VIII/3, 1992; Aagaard et al., 1996; SWERUS Scientific Party, 2016). Sediments were dry-sieved using a ≥150-µm sieve and all samples were randomly split into ∼1/32 of the original volume using a microsplitter before sorting by size fraction. The samples were then sorted into <150-µm, 150–212-µm, 212–350-µm, and >350-µm size fractions. We obtained radiocarbon dates by picking one milligram of N. pachyderma morphotype 2 (Nps-2) from the 150–212-µm fraction, which were measured at the Keck Carbon Cycle Laboratory at the University of California, Irvine. We calibrated radiocarbon ages using CALIB REV8.2 with a marine reservoir correction of 238 years (based on Marine20 data) and, using an uncertainty value of 70 years, calculated from the most proximal Arctic Ocean sites in the database (Stuiver & Reimer, 1993; Heaton et al., 2020). All ages in this study are reported in calibrated years before present (1950). The ages of five core-top samples ranged from 0.3–1.5 ka (SWE13, SWE14, SWE18, SWE22, and B31), one sample was dated to 3.4 ka (PS2170), and three others were between 5–5.8 ka (B16C, B20A, and PS2185). We omitted one core top older than 6 ka (not shown in figures) to constrain the temporal range of this study and provide a mid- to late Holocene baseline of foraminiferal assemblages and geochemistry.

Across our sampling domain, the upper 100 m of the Arctic Ocean is relatively cold and exhibits low spatial variability in mean annual temperatures at 50 m (Fig. 1a; see Fig. A1 for contextualized bathymetric features). However, upper Arctic Ocean waters are comprised of a mixture of Pacific and Atlantic waters, as well as continental freshwater input, creating a steep surficial salinity gradient and halocline (Timmermans & Marshall, 2020). Four of our sites are located in the Eastern Siberian Sea where the Kolmya and Indigirka rivers terminate onto the continental slope. Adjacent to the Eastern Siberian Sea, in the Canadian Arctic (Fig. A1), is the much less saline Beaufort Gyre (Fig. 1b), where freshwater influx, circulation conditions, and the formation of sea ice modify sea surface salinity gradients and stratification throughout the year (Timmermans & Marshall, 2020). Although we do not have samples from the Beaufort Gyre, we collate previously published measurements proximal to the Beaufort Sea that allow us to distinguish between the low- and high-salinity regimes (Spielhagen & Erlenkeuser, 1994; Xiao et al., 2014).

Five other samples mark a transect across the central Arctic. Toward the Fram Strait, more saline North Atlantic waters enter the Arctic, influenced by the North Atlantic Current, and flow deeper than North Pacific waters and freshwater influx in the water column (Cronin et al., 2014; Timmermans & Marshall, 2020). The percent contribution of Pacific and Atlantic water to the Arctic water column is correlated with the relative depth of the halocline; as the contribution of Pacific water decreases, the halocline shoals (Timmermans & Marshall, 2020). Central Arctic bathymetry includes the Lomonosov Ridge (Fig. A1), which provides a physical boundary for incoming Atlantic water and roughly marks where Pacific water becomes negligible. In contrast, sea ice meltwater increases its relative contribution to the upper 300 m in the central Arctic Ocean, particularly the Amundsen Basin (Ekwurzel et al., 2001). Riverine influx is consistent in its contribution between the Chukchi Abyssal Plain and the Amundsen Basin, however its relative contribution increases in the central Arctic Ocean (Ekwurzel et al., 2001).

Arctic Ocean temperature in the upper 200 m ranges from approximately –2–0°C and salinity from 30–35 (Fig. 1, Figs. A2–A4). Toward the Bering Sea, where Pacific water is the primary component of the upper 200 m, the Eastern Siberian Sea is comprised of two distinct sections due to the strong lobe of riverine freshwater from Siberia. East of 160°E, the water column is primarily comprised of Pacific water with an average salinity of 29.7–32.2 and temperature of –0.17–0.58°C (Semiletov, 2005). West of 160°E, the water column is dominated by riverine influx with an average salinity of 22.3–24.5 and temperature of 2.23–2.63°C (Semiletov, 2005). The average coverage of sea ice in the Eastern Siberian Sea is approximately 90–100% in March and 0–60% in September (Cavalieri & Parkinson, 2012). Sea ice coverage increases moving toward the northern coast of Greenland (Fig. 1c) and the central Arctic has an average coverage of 100% in March and 80–95% in September (Cavalieri & Parkinson, 2012). Overall, our sites cover sharp gradients in sea ice coverage and oceanographic properties.

Given that not all our core tops preserve modern sediments, we assessed the potential for post-depositional alteration of planktic foraminiferal calcite across the sites. We investigated ultrastructure preservation via micro-imaging techniques wherein we imaged individual tests using a Hitachi 3400N Scanning Electron Microscope at the University of Arizona LaserChron Laboratory. Samples were uncoated and photos were taken at 200–500×, 750×, and 1500× to image the whole test and the underlying ultrastructure. All images are secondary electron images except those of the unencrusted foraminiferal outliers (UFOs; Fig. 2 and Figs. A9–A11), which were backscattered electron (BSE) images. A total of 23 specimens were imaged across the nine core-top samples, where individual tests were randomly selected from assemblages. We evaluated (N. pachyderma morphotype) Nps-2 images across all sites for relative dissolution using preservation descriptions outlined by Romanova et al. (2017) and Consolaro et al. (2015), shown in Figure A5. Specimens were graded on a scale of 1–5, where a grade of 1 consists of total denudation or removal of surface crystals on test crust, large pores, noticeable thinning of ultimate chamber, channeling of aperture lip, and evidence of recrystallization. A grade of 5 consists of sharp pyramidal surficial crystals, an intact aperture lip, thick-walled ultimate chamber, and small pores. Grades between 1 and 5 represent the relative presence of dissolution or recrystallization characteristics (Fig. A5). Additionally, we examined the benthic-planktic ratio as a proxy for preferential dissolution of planktic species between sites (Thunell, 1976; Fig. A5).

We generated foraminiferal assemblages (average n = 317 specimens) at each of our sites and enumerated the relative planktic foraminiferal species abundances and morphotypical abundances of N. pachyderma (Figs. 3–4). Abundances were calculated from a randomly split sample fraction prepared prior to size-dependent sorting. We prepared assemblages from randomly sorted fractions of bulk sedimentary core material greater than 150 µm. Planktic species and morphotypes of N. pachyderma were identified following Altuna et al. (2018), Schiebel & Hemleben (2017), and Kennett & Srinivasan (1983). Neogloboquadrina pachyderma morphotype assemblages were recorded from the same species assemblages. We labelled specimens as “ambiguous” if those individuals contained hybrid characteristics of two distinct species, kummerforms, or other aspects of ontogeny and morphotypical variability that precluded a clear demarcation of species. For example, ambiguous N. pachyderma individuals included those with secondary individuals attached to the final chamber (Fig. A11), tests with extra half-chambers in the final whorl, and tests with malformed ultimate chambers. Morphotypes of N. pachyderma are labeled Nps-1, Nps-2, Nps-3, Nps-4, and Nps-5 in accordance with Altuna et al. (2018). Following Altuna et al. (2018), we do not employ the distinctions “encrusted” and “nonencrusted” seen in previous work, as ultrastructure characterization is incorporated into morphotypical categorization. Deferment to morphotypical description is also more specified than “encrusted” and “nonencrusted”. Below we provide a brief, colloquial synopsis of morphological traits used in the identification of the four most abundant species in conjunction with the established taxonomy of Altuna et al. (2018), Schiebel & Hemleben (2017), and Kennett & Srinivasan (1983):

• Neogloboquadrina pachyderma (Fig. 2.3–7) typically has four globular chambers of comparable size in the final whorl with a prominent aperture lip. However, there is significant morphotypical expression within the species. Nps-1 has no distinct chambers in the final whorl or the spiral side. It is heavily encrusted with a pit-like umbilical aperture approximately centered on the umbilical side. On some occasions sutures may be visible on less encrusted forms. A half-chamber is sometimes present in the final whorl; it may be described as “egg-like”. Nps-2 is the typical form of N. pachyderma. The form appears quadrate with four chambers in the final whorl and distinct sutures. The aperture has a low arch and is confined to the width of the ultimate chamber. Nps-3 is more rhombic in form than Nps-2 and typically has four- and one-half chambers in the final whorl. The half-chamber is an extension of the ultimate chamber and thus creates an elongated aperture with an irregular arch. Nps-4 has five chambers in the final whorl with a pronounced aperture lip. The aperture is typically more open than the typical form and may appear asymmetric if the last two chambers in the final whorl are fused. Nps-5 is typically rhombic with four chambers in the final whorl and a thinner, but still pronounced, aperture lip. The slight low trochospiral form causes the aperture to open slightly toward the periphery yet remain visible from the umbilical side, resulting in a slightly asymmetrical arch and a less encrusted appearance than in other mature morphotypes.

• Neogloboquadrina incompta (Fig. 2.1) was identified as any test having the same morphological characteristics as N. pachyderma except for its dextral coiling direction (Darling et al., 2006; Schiebel & Hemleben, 2017). We note that we observed large variability in morphology amongst N. incompta individuals as well, but we did not discriminate between them for this study (Fig. A10).

• Turborotalita quinqueloba is a species often cited in marginal Arctic planktic populations (Carstens et al., 1997; Pados & Spielhagen, 2014). Turborotalita quinqueloba typically has five chambers in the final whorl with curved sutures and a slight positive gradation in size toward the penultimate chamber in the final whorl. The ultimate chamber is distinct from the other chambers in the final whorl because it is not globular, but rather ampullate and extended toward the aperture. The lip is pronounced, flattened, and extends over the aperture. The lip sometimes has a wavy edge. This hanging lip creates the appearance of a slit-like pseudo-aperture open toward the edge of the test. Paleo samples include spines. There are many occurrences of individuals which display some of these characteristics along with many characteristics of N. pachyderma. We suggest these may be hybridized individuals. These indeterminate forms are labeled as ambiguous.

• We term a group of planktic specimens as unencrusted foraminiferal outliers (UFOs; Fig. 2.2, Fig. A9), which are characterized by their ovate ultimate chamber and arching aperture open to the edge of the test. Their ultrastructure is atypical for mature neogloboquadrinids in the Arctic Ocean and typically consists of individual conical growths of calcite with better preservation of shape and length near the aperture (Fig. 2.2, Fig. A9). Unlike mature neogloboquadrinids, calcite crystals do not typically grow together to provide uniform test coverage. Under a stereo microscope, the texture appears pearlescent rather than sugary. The final whorl typically has four to five globular chambers with straight sutures. The spiral side has defined chambers that increase in size in a low trochospiral. UFOs have thinner walls than other high latitude specimens and they typically grow to 100–200 μm.

Additionally, rare specimens of Turborotalita humilis, Globigerinoides ruber, and Globigerina bulloides were also identified in the assemblages. To ensure assemblage size did not impact the relative abundance of species, we duplicated three site assemblages (SWE14, PS2185, PS2170) with a difference of at least 250 individuals between replicate and original assemblages (i.e., ∼300 relative to ∼550 specimens). We found no significant differences between relative abundances of duplicated and original assemblages. We thus conclude that our assemblage sizes do not affect the observed relative percent abundance of species or morphotypes.

We measured δ¹⁸O and δ¹³C on Nps-1, Nps-2, Nps-3, Nps-5, N. incompta, and UFO tests, using a range of 30–80 tests from the 150–212-µm fraction for each group (Figs. 4–6). All Nps-2 measurements were duplicated whereas other morphotype and species abundances were too low for replication. We combined analytical uncertainty and the differences between replicates (which were comparable to analytical precision) using root mean square errors. Neogloboquadrina pachyderma-4 was excluded from geochemical analyses due to a lack of available material. For Nps-2, approximately 100 individual Nps-2 tests were picked from the 150–212-µm dry fraction, gently crushed with glass plates, and homogenized using a small paintbrush. Thirty to fifty micrograms of material were partitioned for stable isotope analysis. Stable oxygen (δ¹⁸O) and carbon (δ¹³C) isotope values are reported in δ notation (in ‰) relative to the Vienna Pee Dee Belemnite (VPDB) scale and were measured on an instrumental setup composed of a Thermo Kiel IV Carbonate Device coupled to a Thermo Scientific MAT 253+ Mass Spectrometer housed at the Paleo² Laboratory at the University of Arizona. Based on repeated measurements of IAEA-603, the analytical precision over the sampling interval was 0.07‰ for δ¹⁸O and 0.05‰ for δ¹³C (n = 15), consistent with the long-term precision (1σ of δ¹⁸O = 0.05‰ & 1σ of δ¹³C = 0.03‰) of this setup.

We also collated previously published measurements of δ¹⁸O and δ¹³C on N. pachyderma tests from recent Arctic Ocean sediments (Fig. 6). We sub-selected only those core tops that were dated to be younger than 6 ka (Zahn et al., 1985; Poore et al., 1999; Volkmann & Mensch, 2001; Nørgaard-Pedersen et al., 2003; Hillaire-Marcel et al., 2004; Adler et al., 2009; Polyak et al., 2009; Werner et al., 2013; Xiao et al., 2014; Zehnich et al., 2020). We note that these previous studies did not differentiate between N. pachyderma morphotypes. The average age across the collated samples is 2.32 ka. The average age of the samples analyzed in our study is 2.8 ka.

Forward-modeled foraminiferal calcite δ¹⁸O values were averaged across a one-degree square around core site coordinates prior to comparison with measured values. This combination of equations provided the most reasonable estimates of N. pachyderma δ¹⁸O as compared to other experimental light regimes and species (not shown) presented by Bemis et al. (1998).

We measured Mg/Ca in N. pachyderma tests as a proxy for oceanic temperatures. For elemental measurements, we took the remnant split of the crushed Nps-2 tests and cleaned samples following established Mg/Ca cleaning protocols (Barker et al., 2003). Prior to cleaning, crushed samples were visually inspected to ensure that non-calcareous grains were absent. The cleaning procedure we followed included clay removal using water and methanol rinses, removal of organic matter using a buffered solution of hydrogen peroxide, and ultimately, dissolving for elemental analysis following a weak nitric acid leaching step (Barker et al., 2003). We analyzed samples on a Thermo iCAP 7400 inductively coupled plasma-optical emission spectrometer housed at the Paleo² Laboratory at the University of Arizona. An internal gravimetric standard was employed between samples to correct for instrumental drift and to assess precision (Schrag, 1999). Repeated analysis of external matrix-matched standards using the same analytical grade nitric acid solution yielded a precision of ±0.12% (1σ = 0.018 mmol/mol).

Using the paired δ¹⁸O-Mg/Ca measurements on Nps-2 specimens, we computed temperature and salinity using previously published calibration equations (Fig. 7) and obtained uncertainty estimates using the Paleo-Seawater Uncertainty Solver (PSU Solver) algorithm (Thirumalai et al., 2016). The PSU Solver allows user input of the relationships between foraminiferal δ¹⁸O, Mg/Ca, and temperature and δ¹⁸O_seawater, and iteratively solves for the latter two variables and their respective uncertainties in a Monte Carlo framework (Thirumalai et al., 2016). We reconstruct seawater temperatures using the following published calibration equations, where Mg/Ca is the magnesium to calcium ratio in measured foraminiferal calcite, and T is temperature of calcification:

No Mg/Ca-temperature calibrations have been developed at sites proximal to our sampling region. Furthermore, a range of non-thermal influences including salinity, bottom water calcite saturation, and pH can modulate foraminiferal Mg/Ca and thereby affect calculations of temperature (Khider et al., 2015; Gray & Evans, 2019; Tierney et al., 2019; Holland et al., 2020). However, culture and core top studies have shown that pH and/or carbonate ion concentrations and bottom water calcite saturation do not significantly affect N. pachyderma Mg/Ca (Davis et al., 2017; Tierney et al., 2019). Moreover, recent culture experiments indicate very low N. pachyderma Mg/Ca sensitivity (∼0.5 mmol/mol per >20 salinity changes) to salinity (Bertlich et al., 2021). Equations incorporating the influence of salinity on Mg/Ca, but developed for other planktic species (Kısakürek et al., 2008; Tierney et al., 2015) yielded unrealistic temperatures for our N. pachyderma Mg/Ca values.

Using our paired foraminiferal δ¹⁸O-Mg/Ca data at each site as input into PSU Solver, we calculated δ¹⁸O_seawater (Eqn. 1) and salinity (Eqn. 2) and their uncertainties. We then compared the inverted temperature and salinity estimates to observations from WOA (Fig. 7). We submit that this exercise is qualitative and for illustrative purposes; as the core tops are not all modern, we caution that differences between observations and inverted values may arise due to climatic change from the age of the core top relative to present-day observations.

Additionally, we applied the ‘BayFox’ algorithm to our N. pachyderma δ¹⁸O data. ‘BayFox’ is based on a Bayesian calibration of the δ¹⁸O paleothermometer tailored to reconstruct seawater temperature from stable oxygen isotope data for planktic foraminifera (Malevich et al., 2019). We used Equation 2 to derive seawater δ¹⁸O from WOA observations of salinity from 0–150 m as input into BayFox. We then used 0°C, the approximate temperature of ∼50–100 m depth at our sites, as the prior mean seawater temperature (Timmermans & Marshall, 2020). Finally, we use the standard deviation of WOA seasonal temperature variability at 0–150 m depth at all sites (0.34°C) as the prior on the standard deviation.

Radiocarbon ages confirm that all sediment samples in this study are younger than 6 ka (Fig. 6). We note that our core-top samples are, on average, younger than those used to reconstruct modern conditions using N. pachyderma in previous studies (Fig. 6). Of note is site B31, which is dated to 1.4 ka, exceptionally young for central Arctic Ocean sites. Our radiocarbon dates also indicate that our Eastern Siberian Shelf (ESS) sites (SWERUS 13–22) are significantly younger than our central Arctic core tops as well as Canadian and central Arctic cores from previous studies.

Ratios of benthic to planktic foraminiferal abundances are consistently low across sites (average ≃11%; Fig. A5) and provide confidence for good preservation of planktic individuals. Site SWE22 has an unusually high benthic/planktic ratio (50.7%). The presence of UFOs, a thin-walled group of individuals, further suggests that samples are minimally affected by dissolution. Visual analysis of Nps-2 indicates intact calcite preservation and minimal to no recrystallization (Fig. 2, Fig. A5). We additionally note that the average grading for post-depositional alteration at our site based on microimaging was 3.2, where 5 represents a pristine sample (see Materials and Methods). We note that site SWE18 is notable for low amounts of (thin-walled) UFOs and that this site showed the highest benthic to planktic ratio, indicating the potential for some dissolution. However, the shallow depth (349 m) of this site could have some influence on these observations. We note that other SWERUS sites do not exhibit the same high benthic-planktic ratios and show better preservation and notable percentages of UFOs.

We found no significant correlations between radiocarbon ages and N. pachyderma Mg/Ca or δ¹³C across its morphotypes. Interestingly, we find significant negative relationships between δ¹⁸O across all measured species and the radiocarbon age of our core tops (e.g., Nps-2 depicted in Fig. 6), indicating that younger samples are associated with higher δ¹⁸O values (i.e., consistent with cooler and/or higher δ¹⁸O_seawater/salinity conditions). Dissolution (with or without recrystallization) would serve to preferentially heighten older samples’ δ¹⁸O values (Lohmann, 1995)—something we do not observe. Additionally, we find weak, negative correlations between average δ¹⁸O values and water depth, with shallower sites displaying higher δ¹⁸O. This observation further suggests that our measurements are not affected by dissolution or recrystallization, both of which would induce anomalously positive δ¹⁸O values at greater depths (Edgar et al., 2015).

We find similar relative abundances of planktic foraminiferal species across sites (Fig. 3a). All planktic species percentages are reported relative to the total number of planktic foraminifera. Neogloboquadrina pachyderma is the dominant species, making up at least 85% of the planktic foraminiferal population at each site and 90% on average. Neogloboquadrina incompta (2–7% of all planktic foraminifera, average = 5±0.01%) and UFOs (0–4% of all planktic foraminifera, average = 2±0.02%) are also significantly present at all sites, except for SWE18. Turborotalita quinqueloba is present at every site, however it comprises 0.2–1.4% (average = 0.6±0.00%) of all planktic foraminifera. This contrasts with previous work in the Fram Strait showing a higher prevalence of the species (Pados & Spielhagen, 2014; Schiebel et al., 2017). Benthic foraminifera make up 0–34% of our foraminiferal assemblages with an average of 8±0.12%, relative to the sum of both planktic and benthic foraminifera (Fig. A5). Ambiguous individuals are minimal and make up no more than 4% of the planktic foraminiferal population (Fig. 3).

Morphotypical variability of N. pachyderma was also found to be consistent across sites (Fig. 3b). The most abundant morphotype of N. pachyderma is Nps-2, making up at least 44% (average = 59±0.05%) of N. pachyderma at our sites. The next most abundant are Nps-3 (10–20%, average = 15.5±0.03% of total N. pachyderma) and Nps-1 (7–21%, average = 15±0.05% of total N. pachyderma). Neogloboquadrina pachyderma-4 (2–6%, average = 3.7±0.02% of total N. pachyderma) and Nps-5 (1–9%, average = 4.8±0.02% of total N. pachyderma) are both uncommon, but present. Ambiguous N. pachyderma specimens (see Materials & Methods) make up 0.9–6.1% of the morphotypes of N. pachyderma (average = 3±0.02%). We observe the same morphotypical forms defined in Altuna et al. (2018) in the N. incompta population, however we do not quantify the morphotypical variability.

We find that the average δ¹⁸O and δ¹³C of all analyzed neogloboquadrinids across sites (including four N. pachyderma morphotypes and N. incompta) are similar (Fig. 4). UFO δ¹⁸O values were consistently lower at each site (Fig. 4). Inter-site averages of morphotypes Nps-1 δ¹⁸O (2.35±0.6‰), Nps-2 δ¹⁸O (2.22±0.28‰), Nps-3 δ¹⁸O (2.32±0.39‰), Nps-5 δ¹⁸O (2.07±0.45‰) were statistically indistinguishable within uncertainty relative to each other and relative to average N. incompta δ¹⁸O (2.22±0.45‰). We found no systematic offsets between morphotypes or species within each site, as delineated by the average differences relative to Nps-2 values, the dominant morphotype; differences between Nps-2 and Nps-1 (spatial average = 0.13±0.43‰), Nps-3 (0.11±0.26‰), Nps-5 (–0.15±0.35‰), and N. incompta (0.01±0.34‰) showed no consistent offsets. In contrast, average UFO δ¹⁸O (1.44±0.61‰) differs significantly from N. pachyderma and N. incompta. The difference between Nps-2 and UFO δ¹⁸O (average inter-site difference = –1.1±0.87‰) is consistently negative at each site. Values of δ¹³C for all measured N. pachyderma morphotypes, N. incompta, and UFOs are statistically indistinguishable (Fig. 4).

To obtain insights into average depths of N. pachyderma calcification across our samples, we compared Nps-2 δ¹⁸O with forward-modeled δ¹⁸O at various depths (Fig. 5). We stress this exercise is not an ideal comparison as our core tops are not modern and four of our core tops are older than 1,500 years. Nevertheless, we find that N. pachyderma δ¹⁸O at our youngest sites located in the Eastern Siberian Sea compare well with forward-modeled δ¹⁸O using annual mean temperature and salinity observations at 50–100 m depth (Fig. 5). However, the forward model appears to overestimate measured δ¹⁸O_calcite toward the central Arctic, where samples are, on average, older than 1.5 ka and are apparently consistent with near-surface temperature and salinity.

Our N. pachyderma stable isotope measurements fill in key spatial gaps for Arctic core-top studies (Fig. 6). Although our values are within the range of those from previous studies of mid- to late Holocene sediments (Fig. 6b–c; note that they did not discriminate for morphotypes), we also observe distinct differences imposed by regional Arctic oceanography. Across new and previous measurements, we find that average mid- to late Holocene foraminiferal δ¹⁸O is 1.91±0.52‰ and δ¹³C is 0.95±0.22‰ across the Arctic Ocean. Average N. pachyderma δ¹⁸O is slightly higher in our study (average δ¹⁸O = 2.09±0.24‰) than previous studies (average δ¹⁸O = 1.79±0.30‰) in the central Arctic Ocean (1500–3000 km from the Bering Strait). This slightly lower average is driven by low δ¹⁸O values in the older cores closer to the Bering Strait (Fig. 6b). There is a clear deviation between δ¹⁸O in our sites and previous work in the Pacific-adjacent region of the Arctic Ocean (Fig. 6b). Limiting the distance from the Bering Strait to 0–1600 km, N. pachyderma in the Eastern Siberian Sea are enriched in ¹⁸O (average δ¹⁸O = 2.47±0.18‰) relative to samples from the Chukchi Sea region (average δ¹⁸O = 1.56±0.24‰). We note sites next to each other show similar isotopic values, denoting different isotopic regimes within the data. Average N. pachyderma δ¹³C from this study (0.84±0.12‰) is indistinguishable from δ¹³C of previous studies within observed intersite variability (0.98±0.23‰; Fig. 6). We find that spatial variations in N. pachyderma δ¹³C are much smaller in magnitude across the Arctic Ocean than compared to their δ¹⁸O composition (Fig. 6).

We assessed the suitability of several Mg/Ca-temperature relationships for our core-top measurements of N. pachyderma (Table 4). Jonkers et al. (2013) developed a temperature equation using North Atlantic Ocean sediment traps over a temperature range of 5–10°C. The equations of Kozdon et al. (2009) and Nürnberg et al. (1996) were established using sub-Arctic samples from the Norwegian Sea, with temperature ranges of 3–6°C and 0–15°C, respectively. Vázquez Riveiros et al. (2016) used Southern Ocean samples to develop an equation with a temperature range of –1–9°C. Livsey et al. (2020) utilized samples from Fram Strait plankton tows and cultured specimens captured in Bodega Bay to develop a Mg/Ca-temperature equation within –2–12°C. When compared with other equations, Livsey et al. (2020)’s relationship provides the most reasonable estimates for Arctic sites relative to observed temperatures (Table 4). All other equations yield average inter-site temperatures ranging from 5–11°C—unrealistically warm for Arctic sites, across any part of the water column (Fig. 7).

We then compared the reconstructed temperatures with the temperatures observed at depth. Conductivity and temperature at depth (CTD, Fig. 7) observations were available at locations close to our core sites (Anderson, 2006; Swift, 2006; Rabe & Wisotzki, 2010; Björk, 2017). Across these sites, observed temperatures at 0 m, 50 m, 100 m, 150 m, and 200 m range from –1.80 to 1.09°C (average = –1.12±0.89; Fig. 7). From 0–50 m (of these selected depths) temperatures range from –1.80 to –0.98°C (average = –1.57±0.21). From 0–100 m, temperatures range from –1.49 to 1.90°C (average = –0.16±1.02). From 0–50 m, temperatures are relatively constant from site B31 to site SWE13 (from the central Arctic Ocean to the Eastern Siberian Shelf); however, at 100–200 m, temperatures slightly decrease from B31 to SWE13 and exhibit higher variability (Fig. 7). World Ocean Atlas temperatures at these sites are similar to CTD observations, but slightly underestimate temperature from 0–100 m by an average of 0.03±0.4°C. At 150–200 m, this underestimation is more noticeable (average = –0.28±1.06°C), but still less than 1°C on average. Temperatures calculated from Nps-2 Mg/Ca using the Livsey equation are relatively warmer compared to observed temperatures, which range from –2–0.5°C in the upper 200 m of the Arctic Ocean when excluding the Fram Strait (Pnyushkov and Polyakov, 2022; Fig. 7). Nevertheless, compared to reconstructed temperatures using the Livsey et al. (2020) equation, we find that N. pachyderma Mg/Ca-temperatures at most sites suggest calcification at 150 m or deeper, with the East Siberian Shelf sites exhibiting the warmest reconstructed temperatures and, by extension, deepest reconstructed calcification depth (Fig. 7). At these sites, which are more recently aged (<1.5 ka) than others and likely compare better with observations, the Livsey et al. (2020) equation appears to predict temperatures that are too warm, especially considering that forward-modeled δ¹⁸O predicts shallower (∼50–100 m; Fig. 5, Fig. A7) apparent calcification depths.

Reconstructed δ¹⁸O_seawater values across core tops ranged from –2 to 0‰ (VSMOW). We note that these compare well with proximal δ¹⁸O_seawater measurements from the Laptev Sea, which range between –2.14 and 0.22‰ (VSMOW) at water depths of 40 to 200 m (Bauch & Cherniavskaia, 2018). Salinity calculated from our core-top δ¹⁸O_seawater values (Eqn. 2) ranged from a maximum of 34.8±0.5 at site SWE14 to 31.6±0.7 at site B31 (Fig. 7). We then compared these values to WOA salinity at various depths (Fig. 7b), which showed a consistent trend of deepening toward the East Siberian Sea sites, in agreement with forward-modeled foraminiferal δ¹⁸O (Fig. 5). However, we acknowledge that deeper apparent calcification depths of the East Siberian Sea according to the salinity comparison (∼100–200 m; Fig. 7b) relative to forward-modeled calcite δ¹⁸O (∼50–100 m; Fig. 5) may arise due to uncertainty in calibrations of foraminiferal Mg/Ca composition and seawater temperature. Similar to expectations from reconstructed temperature, this exercise also suggests that N. pachyderma in the Eastern Siberian Sea region may calcify at deeper depths than individuals in the central Arctic (Fig. 7b). We further compared our measurements with seasonal trends in forward-modeled δ¹⁸O, temperature, and salinity (Figs. A7–A8) but found that seasonality amongst sites was unable to explain our observed trends.

Our work provides a mid- to late Holocene baseline of planktic foraminiferal species abundances and benthic-planktic ratios, N. pachyderma morphotypical abundances, δ¹⁸O and δ¹³C of multiple foraminiferal species, and Nps-2 Mg/Ca measurements across nine East Siberian Shelf and central Arctic core tops. The youngest core top dates to 0.37 ka, five core tops are younger than 1.5 ka, and four range from 3.36 to 5.84 ka. The older core-top ages occur due to relatively low sedimentation rates and dynamic sedimentary settings across Arctic subregions (Bröder et al., 2016; Gemery et al., 2017), and are a consistent source of complexity in Arctic core-top studies (Xiao et al., 2014). Bioturbation in the Arctic Ocean is relatively low, with mixing rates ranging from ∼0.27 cm/yr in deep water sites near Svalbard to ∼0.04 cm/yr in the central Arctic Ocean under permanent ice (Soltwedel et al., 2019). The mid-Holocene core tops in our samples preclude quantitative comparisons of measurements with modern observations across all sites. As such, we refrain from attributing precise depths of N. pachyderma calcification or performing Mg/Ca-temperature calibrations and contrast all sites with observations for illustrative purposes (Figs. 5, 7). Despite the older core tops, the lack of correlation between water depth and δ¹⁸O, δ¹³C, and Mg/Ca, respectively; consistently low benthic-planktic ratios across samples; and microimaging analysis collectively point to the pristine nature of calcite preservation in our samples.

All sites have approximately the same distribution of planktic foraminiferal species and N. pachyderma morphotypes. Despite the observed surface and subsurface temperature and salinity gradients across our sites (Figs. 1, 7), we find no obvious relationships between local oceanography and relative distribution of species and N. pachyderma morphotypes. Although we do not observe modern spatial trends in morphotypical distribution, we do not discount N. pachyderma morphotypical assemblages as a tool for studying major changes in Arctic paleoceanography. Even though there are sharp changes in salinity and sea ice distributions across these sites, morphotype and species assemblages may be responding to temperature gradients, which are much smaller than salinity changes, and the average makeup of Pacific versus Atlantic waters at these sites. These oceanographic parameters might have changed significantly over the Quaternary (Poore et al., 1999), and we therefore suggest that spatially distributed downcore studies of assemblages may provide more insight into this issue. Overall, however, our measurements suggest relatively stable planktic foraminiferal species and N. pachyderma morphotype distributions from the Eastern Siberian Shelf to the central Arctic and provide a baseline for future downcore work.

Previous work has discounted the presence of N. incompta in the Arctic and surrounding areas (Altuna et al., 2018; Volkmann, 2000). However, we find that they are significantly present across our sites, exhibit the same morphotypical diversity as N. pachyderma, and exhibit virtually identical stable isotopic values as Nps-2. Neogloboquadrina incompta occurs more frequently than 1–3% of the planktic assemblage, indicating it is likely not a type of N. pachyderma (Darling et al., 2006), and is distinct at the species level; future genetic analyses can confirm or refute this hypothesis. Darling et al. (2006) also notes N. pachyderma is almost entirely sinistral at high latitudes (low temperatures) and shows “type N. incompta” plates which confirm our observations, adding further confidence to our identification. We follow convention and label all dextral neogloboquadrinids adhering to the appropriate taxonomy N. incompta (Darling et al., 2006; Ovechkina et al., 2010; Schiebel & Hemleben, 2017; Lam & Leckie, 2020; Brummer & Kucera, 2022). Furthermore, trends in the relative abundances of N. incompta and N. pachyderma across sites are distinct (Fig. 3). The isotopic similarity of N. incompta to N. pachyderma (Fig. 4) makes it unlikely that these specimens are present solely due to transport from the North Atlantic, as has been asserted previously by Volkmann (2000). We also recognize previous work by Davis et al. (2020) which suggests that high stress environments may cause N. pachyderma to reproduce with a higher percentage of dextrally coiling offspring, although we defer to previous in situ observations from plankton tows supporting the findings of Darling et al. (2006) and the conclusion these are N. incompta rather than N. pachyderma (dextral). Significant changes in abundance and geochemistry of N. incompta should indicate dynamic changes in species’ behavior relative to modern observations.

Similarly, UFOs are consistently present across our sites (Fig. 3), suggesting that the group is not transported from sub-Arctic waters and populates the Arctic Ocean. Additionally, UFO samples exhibit δ¹⁸O values consistent with Arctic oceanographic conditions (Fig. 4), making it likely they are an in-situ population. Unencrusted foraminiferal outliers are the only specimens in this study that consistently show an isotopically distinct regime relative to the other species. Although they are similar to N. pachyderma and N. incompta in δ¹³C, they exhibit relatively lower δ¹⁸O values. Lower δ¹⁸O values suggest a warmer and/or less saline calcification environment, which is found in surface waters at our sites. This suggests that UFOs likely calcify at a shallower depth relative to N. pachyderma and N. incompta. This group is abundant enough to perform geochemical analyses, and we suggest it may be a useful species to reconstruct changes in surface-ocean mixing in the Arctic Ocean.

Unencrusted foraminiferal outliers are likely immature N. pachyderma. The unencrusted texture of solitary, conical calcite crystals shares some similarity with the ultrastructure of Nps-5 in Altuna et al. (2018). However, the geometry of the test is much more planar and rhombic than the quadrate forms of mature N. pachyderma, and the apertures of UFOs are consistently open, ovate, lacking a distinct apertural lip, and open to the edge rather than the umbilical side (Fig. A9). We do not observe clear evidence of spines and suggest their ultrastructure does not resemble sub-Arctic or temperate species. We found several mature individuals in our samples which fit Nps-5 morphology. When comparing stable isotopic values of UFOs and Nps-5, it is clear that Nps-5 values are similar to other neogloboquadrinids, yet UFO values appear as outliers relative to the neogloboquadrinids (Fig. 4), establishing them as geochemically and morphologically distinct. Therefore, we posit UFOs are immature neogloboquadrinids with atypical morphology due to their developmental stage.

This is consistent with the lower δ¹⁸O values observed in UFOs, as we would expect unencrusted, immature individuals or juveniles to calcify higher in the water column and descend as they mature (Hillaire-Marcel et al., 2004). Tell et al. (2022) also supports this understanding of N. pachyderma ecology and vertical migration. We note that we observe UFOs to be a similar size to mature neogloboquadrinids, which aligns with their models of a fixed depth habitat, (i.e., the relative age of the individual does not impact the preferred depth habitat, resulting in a mix of weights and sizes at all depths where individuals are observed. Additionally, it has been shown that N. pachyderma exhibit a large range of morphological variability during asexual reproduction, making it likely immature individuals do not conform to typical mature taxonomy (Davis et al., 2020). However, due to the taxonomic features inconsistent with typical neogloboquadrinid forms, comparable size of these specimens to some mature neogloboquadrinids, and evidence of 4–5 chambers in the final whorl, we cannot rule out the possibility of a hybrid species without further work beyond the scope of this paper. For this reason, we refrain from conclusively identifying this group in this paper.

We found near-identical δ¹⁸O and δ¹³C signatures of N. pachyderma morphotypes across our sites, with no evidence for systematic offsets between morphotypes (Fig. 4). In contrast, across Canadian Archipelago samples, Altuna et al. (2018) reported distinct isotopic signatures for different N. pachyderma morphotypes. They suggested that these differences arose due to different test sizes between morphotypes, which may be representative of different calcification parameters including depth and seasonality. This hypothesis has not yet been tested in the “open” Arctic Ocean. Our dataset suggests that morphotypes of N. pachyderma calcify in similar regimes across Arctic Ocean. We also find no significant differences between stable isotope values of N. pachyderma and N. incompta, indicating a shared and variable Neogloboquadrinid calcification environment across the Arctic Ocean. Therefore, we suggest the potential for bias in N. pachyderma geochemical studies related to morphotypes or spurious N. incompta individuals is minimal. Given the difficulty distinguishing between dextral kummerforms of N. pachyderma and true N. incompta and limited sample sizes in the Arctic Ocean, our results support the geochemical utility of non-selective neogloboquadrinid sampling in downcore reconstructions of Arctic paleoceanography.

Poor preservation likely does not drive the trends we observe in foraminiferal δ¹⁸O. Whereas the negative relationship observed between δ¹⁸O and age could arise due to temporal changes between subperiods of the Holocene or due to spatial differences in average hydrographic conditions across sites, we note that dissolution and recrystallization would induce the opposite trend—older samples would have more positive δ¹⁸O values (Edgar et al., 2015). Previous work using biomarkers and ostracode assemblages have shown cooling starting from 7 ka (Duplessy et al., 2001; Fahl & Stein, 2012; Gemery et al., 2017). While the amount of perennial ice likely increased between the mid-Holocene and present, paleorecords show relatively stable sea ice conditions in the central Arctic since 6 ka analogous to modern observations (Fahl & Stein, 2012; Gemery et al., 2017). Thus, we caution that the overestimation of observations by forward-modeled calcite δ¹⁸O in the central Arctic (Fig. 5) may arise spuriously due to older samples capturing mid-Holocene hydrographic conditions. Yet we note that sample B31, which is located proximal to the North Pole, dates to 1.43 ka and also exhibits consistency with a shallower depth according to forward-modeled δ¹⁸O. Alternatively, if these values are representative of mean-state conditions across the mid- to late Holocene, this trend would still be consistent with deeper calcification in the younger Eastern Siberian Sea samples relative to a shallower N. pachyderma calcification depth in the central Arctic. The more marginal samples we analyze are young and do not pose the same temporal problems. Therefore, we suggest that our geochemical measurements reflect in-situ calcification values that capture mean conditions of the mid- to late Holocene and are minimally impacted by dissolution or post-depositional alteration.

Local oceanographic and bathymetric settings likely drive N. pachyderma δ¹⁸O differences between our measurements and previous studies. This separation is not driven by geological age. In the Chukchi Sea, core tops dated at 2–0 ka are isotopically consistent with older samples (Fig. 6). Comparison of values between the youngest cores across the Chukchi Sea and Eastern Siberian Shelf show the same difference. Eastern Siberian Shelf sites are on average deeper than sites in the Chukchi Sea (Fig. 6). The surface water of the Eastern Siberian Shelf primarily consists of meltwater and river runoff with variable temperatures in summer and high salinity, colder waters in winter due to sea ice formation (Wang et al., 2021). Eastern Siberian Shelf mixed layer water (0–50 m) ranges from –1.7–3.0°C and <28–31 and halocline waters range from –1.7–1.4°C and 29–34 (Wang et al., 2021). The Chukchi Sea mixed layer (0–20 m, 0–6°C and 29–32) and the halocline (–2–0°C and 32–35) are both comparable to the East Siberian Shelf (Clement Kinney et al., 2022). However, the available sites in the Chukchi Sea are closer to the Beaufort Sea and therefore the Beaufort Gyre, which is associated with the multidecadal storage of freshwater sourced from sea ice melt and freshwater influx from North American rivers (Proshutinsky et al., 2019). This creates a strong salinity gradient in the upper water column with more saline conditions in the Siberian Arctic and less saline conditions in the Canadian Arctic. In addition to the oceanographic data, the relatively higher Eastern Siberian Shelf Nps-2 δ¹⁸O suggests a cooler and more saline environment than N. pachyderma in the Chukchi Sea samples. Considering that these comparisons are consistent even among younger sites, we maintain that the ∼1‰ difference between the Chukchi and Eastern Siberian seas accurately reflects spatial differences in regional oceanography.

Our measurements are consistent with previous studies that have demonstrated variable N. pachyderma calcification depths from 50–200 m across Arctic sub-regions (Bauch et al., 1997; Volkmann & Mensch, 2001; Simstich et al., 2003; Greco et al., 2019; Charette et al., 2020; Jonkers et al., 2022; Tell et al., 2022). Comparisons with seasonal observations suggest that N. pachyderma variably calcify over these depths regardless of season (Figs. A7–A8). Surface dwelling N. pachyderma δ¹⁸O from plankton tows in the outer Laptev Sea range from 1–2‰ whereas subsurface N. pachyderma δ¹⁸O range from 2–3.5‰ (Volkmann & Mensch, 2001). Tows from the Nansen Basin additionally show variable N. pachyderma standing stock and abundances across 0–200 m depths, with a preference for shallower calcification habitats in the interior Arctic, away from the Fram Strait (Bauch et al., 1997). Simstich et al. (2003) showed N. pachyderma in the Eastern Greenland Current followed isohalines of 34–35, similar to our observations (Fig. 7). Recently, Greco et al. (2019) hypothesized using plankton tows that sea ice cover influences N. pachyderma depth habitat, with regions of more perennial sea ice associated with shallower calcification habitats and vice-versa (Greco et al., 2019). They further posited that this shoaling underneath sea ice may be related to light availability and average water column productivity (Greco et al., 2019). Previous measurements of N. pachyderma δ¹⁸O in Arctic surface sediments (Fig. 6) also suggest variable calcification depths (Xiao et al., 2014); however, these studies did not measure paired Mg/Ca-δ¹⁸O, which can lead to uncertainty related to the relative influences of δ¹⁸O_seawater and temperature on calcite δ¹⁸O. Our paired δ¹⁸O-Mg/Ca measurements on N. pachyderma provide new insights into N. pachyderma calcification.

At the youngest sites in the Eastern Siberian Shelf (ESS), we find relatively elevated δ¹⁸O (Tables 2, 3; Fig. 6b), Mg/Ca-temperature reconstructions (Table 4; Fig. 7a), and reconstructed salinity (Fig. 7b), consistent with calcification at subsurface depths composed of relatively warm and salty waters. In the ESS, sea ice cover is lower than in the central Arctic (Fig. 1c). Inferred calcification depths at these younger ESS sites differ between δ¹⁸O (50–100 m; Fig. 5), temperature (>200 m; Fig. 7a), and salinity (100–150 m; Fig. 7b). While we recognize modeled δ¹⁸O_calcite at 0 m in the ESS is driven by freshwater input from Siberia, our results suggest it is unlikely that N. pachyderma calcifies within the freshwater-dominated uppermost layers (Figs. 5–7). The development of local Mg/Ca-temperature calibrations will lead to more accurate evaluations of N. pachyderma calcification depths in the ESS. On the other hand, our measurements toward the central Arctic, where sea ice is more perennial than in the ESS, are consistent with a shallower N. pachyderma depth habitat, in line with Greco et al. (2019). Although these samples are older than the ESS samples toward the coast, similar results (and cooler temperatures) at the relatively recent site B31 proximal to the North Pole support a shallower N. pachyderma calcification habitat (0–50 m) in the central Arctic.

In general, marginal Arctic sites exhibit higher Mg/Ca values than at central Arctic locations. We note that the lowest value occurs at site B31 (the geographic North Pole; Table 3). However, two sites do not adhere to this overall trend. SWE22 samples have relatively low Mg/Ca compared to the other SWERUS site samples, and B16C samples display higher values relative to the more central Arctic sites. B16C Mg/Ca could be an anomalous value due to its age, as we would expect anomalously high Mg/Ca in the mid-Holocene relative to modern conditions (McKay et al., 2018). Alternatively, relatively higher Mg/Ca at this site could also arise due to the site’s proximity to sea ice, which creates a much larger range of potential Mg/Ca values for lower temperatures (Vázquez Riveiros et al., 2016). Sea ice melt increases total alkalinity in the surface ocean (Lebrato et al., 2020) which could in turn lead to lower Mg/Ca in both seawater and foraminiferal calcite (Vázquez Riveiros et al., 2016), causing sites with a higher coverage of sea ice (Fig. 1c) to exhibit lower Mg/Ca relative to B16C. Regardless of the exact mechanism, our Mg/Ca values are consistent with measurements (ranging ∼0.5–1.4 mmol/mol) in the Southern Ocean (Vázquez Riveiros et al., 2016). Therefore, it is most likely our measured Mg/Ca accurately reflects seawater conditions at the time of N. pachyderma growth; however, the available (non-local) calibrations cannot accurately capture the empirical relationship between foraminiferal Mg/Ca and seawater temperature in the Arctic.

Neogloboquadrina pachyderma is the dominant species of planktic foraminifera observed in the Arctic Ocean core tops we investigated. Neogloboquadrina incompta was the second most abundant species at our sites. Morphotypes of N. pachyderma and N. incompta appear to calcify in the same local conditions, while unencrusted foraminiferal outliers (which we believe are predominantly immature N. pachyderma individuals) likely calcify in shallower waters. The stable isotopic overlap of N. pachyderma morphotypes points to a similar calcification environment over the mid- to late Holocene. Consistent with previous work, we find, via paired δ¹⁸O-Mg/Ca measurements on N. pachyderma, that this species does not calcify at a singular static depth in the Arctic and instead calcifies at approximately 50–150 m. We find that some Mg/Ca-temperature relationships are more suitable than others for N. pachyderma in our samples and stress the need for more Arctic-specific calibrations. Our data are consistent with the hypothesis that N. pachyderma calcification depth shoals from the continental shelf toward the central Arctic. This pattern coincides with increasing sea ice coverage as well as increasing upper-layer salinity. Relative species and morphotype abundances as well as planktic foraminiferal geochemical measurements presented here should provide a mid- to late Holocene baseline for future investigations of Quaternary climate change in the Arctic Ocean.

We thank Dr. Zachary Michaels for his assistance with imaging samples and acknowledge the help from Dr. Lael Vetter as well as other members of the University of Arizona Paleo² Laboratory. This work was completed with partial funding from the Cushman Foundation; we thank the foundation for their support. We would also like to thank Dr. Jessica Tierney and Dr. Diane Thompson for their feedback and input on this publication. KT acknowledges the University of Arizona Technology and Research Initiative Fund (TRIF) for support, and NSF OCE #1903482, which partially funded this work. Thomas Cronin and Laura Gemery are funded by the U.S. Geological Survey, Climate Research and Development Program and supported by the U.S.G.S. Land Change Program. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The Appendix with Table A1 and Figures A1-A11 can be found linked to the online version of this article.

Table A1. The number of samples relevant to this study per run are reported in conjunction with the number of standards used to calibrate the run and the average stable isotopic ratios of the standards for the run.

Figure A1. An illustration of the primary bathymetric features mentioned in this study relative to core sites. Pink circles represent this study, orange triangles represent previous studies.

Figure A2. a) Average annual sea surface temperature. b) Average annual sea surface salinity. c) Forward-modeled average annual δ¹⁸O_calcite at 0 m.

Figure A3. a) Average annual temperature at 100 m. b) Average annual salinity at 100 m. c) Forward-modeled average annual δ¹⁸O_calcite at 100 m.

Figure A4. Top: a) Average annual temperature at 150 m. b) Average annual salinity at 150 m. Bottom: a) Average annual temperature at 200 m. b) Average annual salinity at 200 m.

Figure A5. SEM images of Nps-2 at each core site. Benthic-planktic ratios and relative grades described in Materials and Methods shown for corresponding images. Sites are arranged by proximity to the Bering Strait.

Figure A6. Output from BayFox. Reconstructed temperature for all sites.

Figure A7. Differences between forward-modeled δ¹⁸O and observed δ¹⁸O from N. pachyderma during a) DJF, b) MAM, c) JJA, and d) SON. Points closer to the bold dashed line at zero indicate more similarity between modeled and measured values.

Figure A8. Seawater temperature (pink circles, top plot, left axis) at each core site (arranged in proximity to the Bering Strait) reconstructed from Nps-2 Mg/Ca data (Eqn. 3). Dashed lines (right axis) show WOA temperature data at 0–200 m for a) DJF, b) MAM, c) JJA, and d) SON. Salinity (blue circles, bottom plot, left axis) calculated using reconstructed temperature (Eqn. 3) and Nps-2 δ¹⁸O (Eqn. 2) for a) DJF, b) MAM, c) JJA, and d) SON. Dashed lines (right axis) represent WOA salinity data. Depths are labeled accordingly.

Figure A9. All backscattered electron images show the umbilical view of unencrusted foraminiferal outliers. Images 1a, 5a, 7a, 12a, and 16a show high magnification views of test ultrastructure. Images 1-3 are from site PS2185. Images 4-15 and 17-19 are from site PS2170. Image 16 is from site SWE22.

Figure A10. Examples of N. incompta morphological variability from sites PS2170 and PS2185. All images utilize back scattered electron methodology.

Figure A11. Examples of kummerforms of N. pachyderma with secondary individuals attached to the penultimate chamber. All images utilize back scattered electron methodology.