Coupled Coccolith-Based Temperature and Productivity High-Resolution Reconstructions in the Eastern Equatorial Pacific During the Last Deglaciation and the Holocene

We present a new high-resolution reconstruction of annual sea-surface temperatures (SSTa) and net primary productivity (NPP) using novel coccolithophore-based models developed for the Eastern Equatorial Pacific (EEP). We combined published coccolithophore census counts from core-tops in the Eastern Pacific with 32 new samples from the Equatorial region, to derive a new statistical model to reconstruct SSTa. Results show that the addition of the new EEP samples improves existing coccolithophore-based SST-calibrations, and allow reconstructing SSTa in the EEP with higher confidence. We also merged the relative abundance of deep-photic species Florisphaera profunda in the same surface sediment samples with existing calibration datasets for tropical regions, to reconstruct annual NPP. Both temperature and productivity calibrations were successfully applied to fossil coccolith data from Ocean Drilling Project Site 1240, in the EEP. The coccolith-based SSTa estimates show a cooling during the Last Glacial Maximum (LGM) and the Younger Dryas, and warming at the start of the Holocene. This pattern differs in the timing and magnitude of the temperature changes from other available SST-reconstructions based on biogeochemical and faunal proxies. We discuss these discrepancies to be the result of different proxy sensitivities to insolation forcing, seasonal bias, and/or preservation artifacts. Reconstructed annual NPP shows a general decreasing trend from the late last glacial period to recent times, which we relate to the weakening of wind-driven equatorial upwelling towards the Holocene. We also calculated carbon export using our SSTa and NPP reconstructions, and compared to other geochemical-based reconstructions for the same location. Our coupled SSTa-NPP reconstruction provides key data to more fully assess the evolution of primary and export productivity as well as organic carbon burial in the EEP, with implications for its role in global biogeochemical cycles across glacial terminations.


INTRODUCTION
The Eastern Equatorial Pacific (EEP) is an important region for Earth's climate. Atmospheric-oceanographic perturbations (i.e. El Niño-Southern Oscillation, ENSO) and extended coastal and upwelling systems occurring in this region influence global climate and biogeochemical cycles in the ocean (e.g., Gu and Philander, 1995;Cane, 1998;Liu and Yang, 2003;Calvo et al., 2011;Costa et al., 2017;Lin and Qian, 2019). The EEP also accounts for 5-10% of the global oceanic net primary productivity (NPP) (Pennington et al., 2006). Mixed layer NPP is inversely correlated to sea surface temperature (SST) in this tropical region (Behrenfeld et al., 2006). When ocean stratification is relaxed, deep mixing brings cooler and nutrient-rich waters to the surface, fueling primary productivity (Behrenfeld et al., 2006;Schneider et al., 2008). The EEP shows large spatial gradients in chemical, physical and biological parameters, subject to strong seasonal and inter-annual variability (e.g., Strutton et al., 2008).
Holocene changes in the NPP and SST dynamics can be inferred from observations and model integrations (Behrenfeld et al., 2006;Timmermann et al., 2014;Krumhardt et al., 2016). These observations show that increasing influence of anthropogenic warming on marine ecosystems has led to a decrease of NPP in tropical regions (Behrenfeld et al., 2006), which is expected to get more pronounced in the coming decades (Laufkötter et al., 2015). Although contentious, recent studies have suggested that perturbations in the oceans could cause a tipping point in tropical SST and NPP in response to global warming, which would have impacts on the interconnected Earth biophysical systems, leading to long-term irreversible changes (e.g., Lenton et al., 2019). However, trends and feedbacks in the climate system are not easy to identify using the relatively short instrumental record, limited to just a few decades of observations (satellites), or a century at most (historical network of in-situ measurements). Palaeoclimatic information derived from sedimentary records provides additional constraints on the response of EEP regional SSTs and NPP to large global climate variations, such as the transition from the Last Glacial Maximum (LGM) to the Holocene (e.g., Leduc et al., 2010;Timmermann et al., 2014;Costa et al., 2017).
Reconstructions of LGM to Holocene EEP productivity also vary between proxies (e.g., Robinson et al., 2009;Schneider et al., 2010;Calvo et al., 2011;Kienast et al., 2013;Patarroyo and Martıńez, 2015;Costa et al., 2017;Diz et al., 2018;Jacobel et al., 2020), with some proxy systems responding to multiple environmental factors (Jacobel et al., 2020;Studer et al., 2021). Sediment cores separated by just a few degrees of latitude, north and south from the Equator, show diverging trends in productivity-related proxies, most likely as a consequence of the strong gradients in nutrient regimes and variable oceanographic features (e.g., Dubois et al., 2010;Dubois and Kienast, 2011;Costa et al., 2017;Studer et al., 2021). Multi-proxy studies thus provide a fuller, and hopefully more accurate, reconstruction of the coupling between stratification and productivity through time, as well as the basis to consider the spatial variation in productivity regimes across the EEP.
Satellite-derived NPP observations calibrated to the relative abundance of the tropical taxa Florisphaera profunda in surface sediments from the Indian Ocean and the South China Sea (Beaufort et al., 1997;Beaufort et al., 2001;Zhang et al., 2016) have been used to quantitatively reconstruct NPP in various regions (e.g., Tangunan et al., 2020). More recently, Hernańdez-Almeida et al. (2019) developed a global calibration comprising samples from all tropical regions, which considers the response of F. profunda to changes in NPP, and allows the reconstruction of NPP across all low-latitude regions with high confidence (e.g., Tangunan et al., 2021).
In this study, we extend these approaches by analyzing coccolith assemblages from a new set of surface sediment samples in the EEP, and merge this data with the dataset from Saavedra-Pellitero et al. (2010), to generate new coccolith-based calibrations to reconstruct both SST and NPP in the EEP region. We then apply both calibrations to the fossil coccolith assemblage at ODP Site 1240, located in the Panama Basin, to provide reliable high resolution NPP and SST quantitative reconstructions for the last 20 kyr covering the late last glacial period, Termination I, and the Holocene.

STUDY AREA AND OCEANOGRAPHIC SETTINGS
The newly compiled dataset of sediment surface coccolithophore assemblages covers an area from 15.7°N to 50.7°S and from 70.5°W to 110.6°W, ranging from the Equator to high-latitudes in the Southern Ocean (see Figure 1), including most of the Peru-Chile Current (PCC) system. The cold, low-salinity water of the PCC flows northward, and the southeast trade winds maintain the PCC and promote the Peruvian and the equatorial upwelling, which makes the EEP a highly productive area ( Figure 1B) (Toggweiler et al., 1991). The PCC waters feed the South Equatorial Current, which flows west driven by the prevailing winds in the region. The North Equatorial Countercurrent represents the most important equatorial surface current flowing east (Tomczak and Godfrey, 1994), which is deflected to the north feeding the North Equatorial Current. The Equatorial Undercurrent flows eastward below the South Equatorial Current and provides nutrient-rich waters to the equatorial upwelling, controlling primary productivity in the region (Murray et al., 1994). The Equatorial Undercurrent continues eastward, feeding the Peruvian coastal upwelling (Toggweiler et al., 1991) and the Peru Undercurrent (Wyrtki, 1981;Strub et al., 1998).
The geographic area that covers the whole surface sediment sample dataset is characterized by a strong gradient in annual mean temperatures, with SST ranging from ca. 28°C north of the Equator, to ca.~20°C in equatorial and Peru-Chile upwelling systems, down to 8°C in southern high latitudes ( Figure 1A). In the EEP, the seasonal changes in SST are linked to shifts in the position of the Intertropical Convergence Zone (Wyrtki, 1966). Upwelling along the Equator and in the eastern boundary current, offshore Peru and Chile, supplies nutrient-rich waters (Berger et al., 1987), resulting in high NPP values (up to 1000 mg C m -2 day -1 ) ( Figure 1B). In contrast, productivity drastically decreases towards the Southern Pacific Gyre where nutrients are very limited. Productivity is lower during boreal winter at the Equator, although winds blowing from the Central American Cordillera cause localized upwelling and decrease of temperatures north of the Equator (Xie et al., 2005).

Modern Coccolithophore and Environmental Data Sets
For the present study we considered new and published data ( Table S1; Supplementary Material). The new dataset consists of 48 surface sediment samples retrieved between ca. 11°N and 12°S (core locations and sample codes can be found in Table S1). The uppermost centimetre was sampled from sediment cores retrieved between 152 and 5086 m water depth. Most of the samples are located above the present Carbonate Compensation Depth -CCD-(4000 m) and Carbonate Lysocline (3800 m) in the equatorial Pacific (Thompson, 1976;Farrell and Prell, 1989). Thirty two surface sediment samples were considered for this study (blue dots in Figures 1A, B), all above the CCD. A further 16 were found to be devoid of coccoliths and therefore excluded from the considered dataset. Based on published radiocarbon dating of planktonic forminifera shells from surface sediments of the EEP (e.g., Koutavas and Lynch-Stieglitz, 2003;Mekik, 2014;Studer et al., 2021) (Table S2; Supplementary Material) as well as on available literature focusing on the south-eastern Pacific (SEP) area (e.g., Abrantes et al., 2007;Saavedra-Pellitero et al., 2010;Saavedra-Pellitero et al., 2011;Saavedra-Pellitero et al., 2013) the age of the sediments can be regarded as no older than late Holocene to Recent, as it mostly lies between 0.79 and 3.57 Cal age ka BP.
The present study follows the coccolithophore taxonomy of Jordan and Kleijne (1994); Jordan et al. (2004), and the electronic guide to the biodiversity and taxonomy of coccolithophores Nannotax 3 (http://ina.tmsoc.org/Nannotax3/index.html) by Young et al. (2020), with few additional considerations, for instance small-sized Gephyrocapsa (< 3 μm) is referred here as small Gephyrocapsa, for consistency with Saavedra-Pellitero et al. (2010). Annual averaged present-day environmental variables such as the annual sea-surface temperature (SSTa) , annual-sea surface salinity (SSS) , apparent oxygen utilization, as well as nitrate and phosphate contents obtained at 10 m water depth in the eastern Pacific Ocean (Garcia et al., 2014) were taken from the World Ocean Atlas 2013 (WOA 13), using Ocean Data View 5.3 (ODV) software (Schlitzer, 2021). NPP estimates were based on MODIS chlorophyll-a, SSTa and photosynthetic active radiation satellite data, using the standard vertically Generalized Production Model (VGPM) (Behrenfeld and Falkowski, 1997).
These environmental data were interpolated to the geographic location of new surface sediment samples, and the ones from Saavedra-Pellitero (2010) to keep consistency in the data, using a weighted average gridding as implemented in ODV. Water depth and distance to the coast were also considered here. The first one refers to the depth at which the surface sediment sample was retrieved, in order to assess potential dissolution effects. The distance to the coast (in km) may control the amount of terrigenous material (including reworked material from continental shelves) delivered to the ocean.

Downcore Dataset: ODP Site 1240
ODP Site 1240 (Leg 202) is located in the southern Panama Basin (0°1.311'N, 86°27.258'W) at 2921 m water depth ( Figure 1B). Studied sediments consist mostly of nannofossil ooze rich in foraminifera with varying amounts of diatoms, and rare siliciclastic components (Mix et al., 2003). The age model for the sedimentary sequence at ODP Site 1240 has previously been established by Pena et al. (2008), and it is based on radiocarbon ages of 11 samples of the planktonic foraminifera Neogloboquadrina dutertrei for the studied interval in the  . (B) Annual mean of net primary productivity (NPP, in mg C m -2 day -1 ) estimated using the VGPM algorithm (Behrenfeld and Falkowski, 1997). The red (Saavedra-Pellitero et al., 2010) and blue (present study) dots indicate location of the surface sediment samples with coccolithophore assemblage data. ODP Site 1240, ME0005-24JC (24JC in the map), MV1014-02-17JC (17JC in the map) and the major surface and subsurface currents in the eastern Pacific Ocean (Strub et al., 1998;Feldberg and Mix, 2002;Fiedler and Talley, 2006)  present work. ODP Site 1240 is ideally located to capture climate and oceanographic signals related to the past upwelling and biological production changes in the EEP, due to the relatively high sedimentation rates (allowing high temporal resolution studies), calcium carbonate content and general good preservation of carbonate microfossils (Mix et al., 2003). Additionally, the availability of multiple climate proxies at the ODP Site 1240 location (or nearby) is crucial in order to improve the reliability of climate reconstructions. The mean sedimentation rate varies from~25 cm kyr -1 during the Termination I to~7 cm kyr -1 during the late Holocene. For this study, a total of 270 samples were collected every centimetre covering the last 19.8 kyr ( Table S3; Supplementary Material).
The temporal average resolution of this study is 66 yr.

Sample Preparation and Microscope Observations
Slides were prepared following the well-established settling technique of Flores and Sierro (1997). Observations and coccolith counts were made with a Leica light microscope at 1000x magnification. A minimum of 500 coccoliths were counted per sample, ensuring that all species with an abundance above 1% were represented (Dennison and Hay, 1967;Fatela and Taborda, 2002). From the 270 samples considered for this study, 127 samples covering the Holocene were already available and therefore included (Cabarcos et al., 2014) (Table S1; Supplementary Material). The coccolith assemblage composition is reported in this work as relative abundance. Additionally, total absolute abundance of coccoliths (N, coccoliths per g of sediment) was calculated with the formula: where n is the number of coccoliths counted in a random visual field; R is the radius of the Petri dish used; V is the volume of buffered water added to the dry sediment; r is the radius of the microscope visual field; g is the dry sediment weight; and v is the volume of mixture pipetted. Knowing the dry density of the sediment and the sedimentation rate, it was possible to calculate the total coccolith accumulation rate (CAR) as follows: where N is the absolute abundance of coccolithophores; d is the dry bulk density of the sediment (https://web.iodp.tamu.edu/), and s is the linear sedimentation rate.

Coccolith Preservation
A dissolution index was calculated to assess the potential effect of preferential carbonate dissolution on the coccolith assemblage. The CEX' index (Boeckel and Baumann, 2004) represents the ratio between small and delicate E. huxleyi + small Gephyrocapsa (= small placoliths) and the strongly calcified Calcidiscus leptoporus as follows: CEX′ = % small placoliths % small placoliths + % C: leptoporus ½Equation 3 CEX' ranges from 0 to 1, with values closer to 1 suggesting a better preservation of coccoliths.

Statistical Analyses
A principal component analysis (PCA) was performed to explore the main gradients in the environmental dataset. The purpose of this analysis was to summarize the multivariate environmental data using components that explain the variance in the dataset (Davis, 1986;Harper, 1999). The 'broken stick' model was used to assess the significance of the PCA axes (Jackson, 1993).
Detrended correspondence analysis (DCA) with detrending by segments (Hill and Gauch, 1980) on the species data was used to estimate the length of the ecological gradient, to decide whether lineal or unimodal multivariate statistical methods should be used (Birks et al., 1990). As the length of the first DCA was shorter than 2 Standard Deviation (SD), i.e. 1.25 SD, we assumed that the species response was lineal, as it was in Saavedra-Pellitero et al., (2013), and applied the appropriate linear ordination method (Ter Braak, 1987;Ter Braak and Prentice, 1988). Redundancy Analysis (RDA) was performed on square root species abundances, with down-weight of rare species and scaling focused on inter-species distances, to explore the environmental-species relationships. We calculated the ratio of the first constrained eigenvalue to the first unconstrained eigenvalue (i.e. RDA1/PCA1 = l1/l2) to assess the significance of each environmental variable on the modern coccolith dataset (Juggins, 2013), with values >1 indicating a strong response of the assemblage to the variable of interest.
For the development of the transfer function, we used the well-stablished regression method Imbrie and Kipp factor analysis (IKFA) (Imbrie and Kipp, 1971), because it is suitable for short gradients. The coefficient of determination (R 2 ) and root mean square of prediction (RMSEP) were assessed by bootstrapping cross-validation (999 permutation cycles) (Vaughan and Ormerod, 2005). The optimal number of factors in the final IKFA model was chosen by establishing a ≥5% improvement in RMSEP over alternatives with less components (Birks, 1998). Outliers were identified as those samples with absolute residuals (difference between observed and predicted values) higher than 2SD of the environmental variable of interest. Similarity between the modern and the fossil dataset was evaluated calculating the chord distance between both (or analogue distance) (Overpeck et al., 1985) ( Figure S1A; Supplementary Material). To evaluate whether the variable in the final model was also significant for the downcore assemblage and ODP Site 1240, we used the statistical test developed by Telford and Birks (2011). This method compares the variance in the fossil assemblage explained by the environmental variable of interest with the proportion of the variance explained by randomly generated environmental variables ( Figure S1B; Supplementary Material). The multivariate ordination analyses, transfer function model, significance and analogue test were performed with the R software (version 3.5; R Core Team, 2018), using the packages 'rioja' (Juggins, 2017), 'vegan' (Oksanen et al., 2019), ggpaleo (Telford, 2020), 'palaeoSig' (Telford, 2019) and 'analogue' (Simpson and Oksanen, 2020). Species thermal niche was calculated adapting the code used in Jonkers and Kucera (2019) (https://github.com/lukasjonkers/species_selection).

Past Primary Productivity Estimates
The percentage abundance of the deep-dwelling species Florisphaera profunda is a well-established proxy for primary productivity (Molfino and McIntyre, 1990b;Molfino and McIntyre, 1990a;Beaufort et al., 1997). The abundance of this species has been calibrated to satellite-derived NPP across multiple oceanic regions, at mid-to-low latitudes (Hernańdez-Almeida et al., 2019) resulting in the relationship: In this calibration core-top dataset the EEP was only represented by 10 samples, taken from Saavedra-Pellitero et al. (2010). Here we re-evaluate and update the F. profunda-NPP relationship in the EEP with the addition of 32 new calibration points. This new calibration is then used to reconstruct NPP for the last 20 kry at ODP Site 1240.

Export Production Estimates
We calculated the export production (e f ) at ODP Site 1240 using the following equation from Laws et al. (2011) and our SSTa and NPP estimates: Where opal is the Th-normalized opal flux and CaCO 3 is the Thnormalized carbonate flux (Kienast et al., 2007).

Distribution of Coccoliths in Surface Sediment Samples From the EEP
Here we focus on core top assemblage data for the EEP sector, comprising the new dataset and a few surface sediment samples located between 15.7°N and 12°S published by Saavedra-Pellitero et al. (2010). There is a wide range of coccolith preservation state recorded across the study area (ranging from poor to good) but for most samples it is good, with a CEX' generally >0.9 ( Figure 2). The best-preserved coccoliths are recovered from sediments in the Panama Basin, whilst the most poorly preserved ones (or samples devoid of coccoliths) are found between~6°N and~12°N,~7°S and~13°S in the EEP region. Gephyrocapsids are the most abundant coccolithophore taxa in the surface sediment samples from the EEP, with small Gephyrocapsa dominating assemblages (35.7% on average) and reaching highest relative abundances > 80% in the southern area of the Panama Basin, at~1°S ( Figure 2). The abundance of small Gephyrocapsa decreases from the Panama Basin towards northern latitudes. Gephyrocapsa oceanica is represented with a mean relative abundance of 23.9%, reaching a maximum of 57.7% offshore Mexico at~15.6°S (Figure 2).
The lower photic zone taxa F. profunda is relatively common in the EEP region, with percentages ranging from 3.6% to 59.7%, and a mean of 26% ( Figure 2). Gephyrocapsa muellerae occurs with a mean relative abundance of 2.6% and maximum of 12.5%, with low relative abundances in the EEP, which are similar to those of E. huxleyi with a mean abundance of 2.1% and a maximum of 14.2% at~0.5°S (Figure 2). The abundance of Calcidiscus leptoporus varies from 0.7% to 9.9%, with a mean relative abundance of 3%. Umbilicosphaera spp. is present with an average relative abundance of 3.4%, reaching maxima of 14.1% at~6.1°S in the Panama Basin, similar as Oolithotus sp. (maxima of 10.9%) ( Figure 2). Other taxa, such as Coccolithus pelagicus, are absent from the EEP, while Rhabdosphaera clavigera is only occasionally present in very low abundances.

Modern Environmental Characterization and Statistical Analyses for the EEP and SEP
The new data generated in this work notably expands the environmental gradients presented by Saavedra Pellitero et al. (2010;2013). For instance, the total SST range considered by Saavedra-Pellitero et al. (2013) was only 6°C, while in the current study is 19°C. The PCA indicates that 65% of the variance is explained by the two main components (PCA1 and PCA2, Figure 3). SSTa and NPP are the main two environmental factors contributing to the PCA1, which explains 35% of the total variance. The PCA shows that the largest environmental gradient corresponds to SSTa, which is strongly anticorrelated to NPP. The rest of the environmental variables initially considered (i.e. SSS, depth, distance to the coast as well as nitrate, phosphate and apparent oxygen utilization) mostly contribute to PCA2, accounting for 25% of the variance.
The RDA shows that there is large compositional difference between the coccolith assemblages from the new dataset and from Saavedra-Pellitero et al. (2010) (Figure 4). However, despite sampling across larger environmental gradients, the species-environment response remains linear (DCA1: 1.25 SD). The variance explained by the environmental variables included in the RDA is ca. 36%. The samples analysed in this study are located in the left side of the ordination plot, related to higher SSTa. Species such as G. oceanica, Oolithotus sp., Umbilicosphaera spp., and F. profunda appear related to these higher SSTa, whilst G. muellerae, E. huxleyi, C. leptoporus, and C. pelagicus appear on the right side of the RDA plot, characterized by lower SSTa and higher NPP values (Figure 4). The ratio between the first constrained axis (RDA1) and the first unconstrained axis (PC1) is the highest for SSTa (1.5), indicating that this variable should be the one used in the model. SSTa explains more than 15% of the coccolith speciesenvironment variation.
The lineal species-environment response and the thermal tolerances of the coccolithophore species indicate that some of the taxa, such as small Gephyrocapsa, show similar abundances across a SSTa range of almost~20°C ( Figure 5). Species such as E. huxleyi, G. muellerae, C. leptoporus and Helicosphaera carteri show higher abundances at SSTa below 20°C, although they do not indicate a uniform increase with decreasing temperatures. Other species show a clearer response to changes in SSTa, including increased abundances of Florisphaera profunda above 14°C, and of C. pelagicus below 13°C. Other taxa show potential thermal niches in the SEP and EEP, ranging 12-17°C for Pontosphaera spp., and 17-20°C for Calciosolenia sp. However, the taxa which show more defined SSTs are not as abundant as taxa with higher ecological plasticity, such as those from the genus Gephyrocapsa.

Coccolith-Based SSTa Model
For the coccolith-temperature based regression, we applied the R 2 of the final model using IKFA regression with 2-factors (assessed using a screeplot) -0.82 -with an RMSEP of 2.4°C ( Figure 6A). There are only minor gaps in the SSTa gradient (9-10°C and 20-21°C) due to absence of surface sediment samples covering these SSTa values in the modern calibration dataset. Interestingly, the fitting line for the EEP (SSTa range between 28-20°C) samples from Saavedra-Pellitero et al. (2010) shows a different slope than the 1:1 line, while the predicted temperatures based on the new EEP surface samples in this study show a better fit to the observed values. No outliers, as defined as residuals higher than 2 * SD the SSTa at 10 m, were identified ( Figure 6B).

Downcore Coccolith Assemblage, Abundance and Preservation
Coccolith CEX' values at ODP Site 1240 are in general above 0.9 suggesting moderate to good preservation ( Figure 8J). CEX' values are below 0.9 only at the beginning of Termination I, between~19 and~17.5 kyr, due to an increase in the relative abundance of C. leptoporus probably due to low SSTs, and perhaps to preservation (not obvious in light microscope observations). The sub-fossil coccolithophore taxa preserved in the surface sediment samples are also present in the section studied. The downcore assemblage is dominated by small Gephyrocapsa, G. oceanica, and F. profunda. Small Gephyrocapsa ranges between 18.4% and 52.5% (mean of 37.5%) and is present in low abundance at the beginning of Termination I, with a marked increase between 17.5 to 16.5 kyr.
High percentages of small Gephyrocapsa are recorded between 16.5 to 15.5 kyr, around 14.5 kyr, and between 13.5 and 12 kyr, followed by a relatively gradual decrease after 13 ka. Gephyrocapsa oceanica makes up a mean relative abundance of 26.4%, and varies from 10.3% to 50.9%. This species is more abundant at the beginning of Termination I with marked drops between 17.5 -17 kyr and during the warm Bølling-Allerød (BA) period, from 14.7 to 12.9 kyr. Lower percentages (ca. 20%) of G. oceanica are observed during the Holocene, with an occasional increase at around 5.5 ka. Florisphaera profunda is an abundant taxon which varies from 7.7% to 45%, and shows increasing values from the beginning of the studied interval towards present times. The coccolith assemblage is dominated by F. profunda during most of the Holocene.
Other taxa are subordinate, with low relative abundances of G. muellerae recorded during Termination I (near 3%), and marked increases around 16.5 and 14 kyr, although relative abundances remained below 1% for many intervals of the studied period. Emiliania huxleyi fluctuates throughout the sedimentary sequence, reaching a maximum abundance of 4.3% around 2.5 ka, and lowest abundance between 11.5 -5 kyr. The relative abundance of C. leptoporus is highest at the beginning of Termination I, especially from ca. 20 -17 kyr (reaching a maximum of 6.9%), and during the late Holocene. Umbilicosphaera spp. reaches relatively high relative abundance before the BA and during the Holocene, up to 5.3%, and has low values between 17 -11 kyr. The total CAR showed maximum values during Termination I with a marked decrease at the end of the YD ( Figure 10A). Increased CAR values were recorded from ca. 10.5 to 5 kyr, and the lowest values during the Holocene, from ca 4.5 to 1.9 kyr.

Temperature and Productivity Reconstructions
We applied the new IKFA-2 calibration model to the downcore coccolithophore assemblages at ODP Site 1240 to reconstruct SSTa (Figure 9). Estimated SSTa vary between~20.7 and 25.4°C for the last~20 kyr, with the lowest values recorded at~19.8 kyr, during the LGM, and highest at the beginning of the Holocene, at~10.5 ka ( Figure 9C). The SST change across the Termination goes from 22.4 to 25.4°C at the highest warming, which is above the calibration error (2.4°C). Other thermal oscillations of smaller magnitude have to be taken with caution, because they are within the magnitude error. A drop in SSTa coincident with the Antarctic Cold reversal, an Antarctic cooling event which started at 14.5 and lasted ca. two millennia (Blunier et al., 1997) is recorded during Termination I, from~14.7 to~11.4 kyr. From the beginning of the Holocene towards modern times (i.e., 1.9 kyr), SSTa fluctuate, but generally decrease (Figure 9). The squared chord distance of the nearest analogue from the modern dataset in this study to each fossil sample for Site 1240 is shown in Supplementary Figure S1A, and "good" analogues could be identified among the modern samples for the fossil assemblages. The analyses for the significance of the environmental variables in the downcore data shows that the fossil assemblage at Site 1240 is driven by SSTa (p<0.05) (Supplementary Figure S1B).
Estimated palaeoproductivity based on the global calibration of the relative abundance of F. profunda (Hernańdez-Almeida et al., 2019) varies between~448.2 and 1071.6 mg C m -2 day -1 ( ± 1.02 mg C m -2 day -1 ) for the last~20 kyr. Highest NPP values are recorded prior to 17 kyr, while low NPP occurs during mid to late Holocene (Figure 10). The reconstructed NPP shows variations with a clear gradual decreasing trend for the last~20 kyr (Figure 10), punctuated by two troughs at 17 and 12 ka.

DISCUSSION
Strong and dynamic changes in temperature and primary productivity occur across the EEP in response to recent major climate Pleistocene-Holocene transitions (e.g., Dubois et al., 2009;Kienast et al., 2013;Costa et al., 2017;Ford et al., 2018). However, conflicting SST and NPP reconstructions through Termination I and the Holocene may reflect real variability in the EEP or might be artifacts of the different proxies and techniques used in various studies (Costa et al., 2017). Almost all proxies are susceptible to seasonal or water depth habitat biases (e.g., Mix, 2006;Dubois et al., 2009;Leduc et al., 2010;Timmermann et al., 2014), but considering these as part of a multiproxy study, alongside the coccolithophore-based approaches presented here, allows a more nuanced assessment of past productivity and temperature variability. Our new models for SSTa and NPP are used to generate a high-resolution reconstruction of oceanographic conditions at ODP Site 1240 over the last 20 kyr.
Assessment of our new IKFA-2 coccolith-based SSTa calibration model shows a better performance (R 2 = 0.82, RMSEP = 2.4) compared to the previous coccolithophore-based transfer   Figure 6A), systematically underestimating SSTs. This trend could be due to the fact that the old approach of Saavedra-Pellitero et al. (2011) was based on Principal Component Analysis and the main factor dominating in the EEP was not considered for the palaeoenvironmental reconstruction. That factor (Factor 1) explained 31.58% of the variance, was correlated with silicate content and SST and dominated by the species F. profunda and G. oceanica. Most likely, the temperature underestimates could be due to the use of G. oceanica in the model, which shows similar abundances across a SSTa range of almost 20°C in the SEP and EEP ( Figure 5). Regarding the follow-up model of Saavedra-Pellitero et al. (2013), low latitude surface sediment samples were excluded in that calibration, therefore the SST range was narrower than in this work and the potential for reconstructions became reduced and restricted just to the SEP. With the addition of the surface sediment samples analyzed in this study, the slope of the predicted versus observed SSTa new model has a value ca. 1 ( Figure 6A), which leads to a more robust calibration in the EEP and allows the model's application to fossil coccolithophore records in the Eastern Pacific, from the Equator to Southern midlatitudes (ca. 47˚S), with increased confidence.

Termination I
All planktonic foraminiferal-based proxies at ODP Site 1240 show a trend of more-or-less consistent warming through the last Termination of~4°C (Lea et al., 2006;Pena et al., 2008;Yu et al., 2012) (Figures 9B, D, E, G). This pattern of warming contrasts with coccolithophore-derived, alkenone-based SST (UK' 37 ) estimates from core ME0005A-24JC (0°01′N, 86°27′W, 2941 m depth) (Kienast et al., 2006), that show a distinct~1.5°C cooling between~18 and 15 kyr, before warming a similar amount by the BA/Younger Dryas (YD) ( Figure 9F). Interestingly, our new assemblage-based SST estimates ( Figure 9C) show an intermediate behavior, with cycles of warming and cooling of magnitudes <1.5°C through the Termination but with only a small long-term background warming trend across the whole Termination. Particular cool intervals are associated with Heinrich 1 (H1) and the YD and, in this, most closely resemble the planktonic foraminiferal Mg/Ca records from the nearby Site TR-163-22 (0°3′N, 92°23′W, 2830 m water depth, Figure 9D) (Lea et al., 2006). We speculate that these patterns could suggest certain degree of seasonal or inter-annual offset between the bloom-forming alkenone-producers, mainly Gephyrocapsa species in this instance, and the wider and phylogenetically more diverse community that contributes to the warmer assemblage-based SSTs (such as Oolithotus sp., Umbilicosphaera spp. and F. profunda; see Figure 5). Additionally, Kienast et al. (2006) suggest that the ENSO-related dynamics could have an overprint on the alkenone-based SST record in the EEP. The offsets between alkenone and assemblage-based SSTs are most marked in the interval between H1 and the Older Dryas (15 -16 kyrs), where alkenone-production appears to be in waters~1.5°C cooler than the temperatures indicated by the whole assemblage.
Increasing seasonality of production through the Termination could also explain the increasing offset between foraminiferal-and coccolithophore-derived SSTs, if algal production was biased towards early upwelling of cold-tongue waters as have been ascribed in the EEP and elsewhere (e.g., Kienast et al., 2006;Leduc et al., 2010;Schneider et al., 2010;Timmermann et al., 2014). For example, foraminifera-based Mg/Ca-data potentially reflects boreal summer insolation changes, as shown by presentday G. ruber fluxes in the Panama Basin (Thunell et al., 1983;Leduc et al., 2010), whilst coccolithophore-derived alkenone SSTs reflecting a bias towards boreal winter temperatures . Taken together, this data indicates seasonal warming of boreal winters through Termination I, but with the persistent presence of colder upwelling waters, with a boreal winter or even Southern Hemisphere-sourced temperature signal (Martıńez-Botı́et al., 2015;Diz et al., 2018).

Holocene SSTs
Foraminiferal Mg/Ca temperatures from both ODP Site 1240 and core TR-163-22 reach peak temperatures in the early Holocene (~6 to 10 kyr), before either cooling slightly and stabilizing (TR-163-22, Figure 9D) or cooling to almost the Recent (ODP 1240) ( Figure 9E). The discrepancies between both Mg/Ca SST records from near-by locations are attributed to the differential impacts of dissolution between the sites (e.g., Mekik et al., 2007;Dubois et al., 2009;Leduc et al., 2010), with foraminiferal preservation getting worse during late Holocene at core TR-163-22 (Lea et al., 2006). In contrast, alkenone-based SST records warm by almost 2°C through the Holocene ( Figure 9F), but only so as to converge on the warmer (~24 - 25°C) Mg/Ca SSTs in the late Holocene. Foraminiferal assemblage SSTs ( Figure 9G) are also consistently warm through the Holocene, but are~2°C warmer than the alkenone SSTs. Again, these patterns suggests a warm-season/boreal signal recorded within the G. ruber Mg/Ca data, whilst alkenone (and foraminiferal assemblage data) track the temperature trends of Southern Ocean-derived upwelling waters. The coccolithophore assemblage-based signal ( Figure 9C) thereby tracks the Mg/Ca data more closely, with a general slight cooling through the Holocene, but with the coolest late-Holocene temperatures of any of the proxies. These trends indicate a long-term warmseason (G. ruber-like) behavior but with a potential cool-bias in the modern calibration data, as indicated by latest Holocene/preindustrial SST estimates of~22 -23°C against measured mean annual temperatures at this location of~24.5°C.

Productivity Reconstruction at ODP Site 1240
In the same way that SST records can vary according to the different proxy sensitivities (Figure 9), conflicting EEP productivity reconstructions for the last 20 kyr may reflect real variability or be artifacts of the proxy techniques used (Costa et al., 2017). Combining a range of proxies responsive to productivity provides a better understanding of the interaction between nutrient availability, consumption and NPP in the surface ocean, its transport in the water column and burial in the sea-floor.

Termination I
While alkenone concentration (Calvo et al., 2011) at ODP Site 1240 ( Figure 10C), show high values during the LGM and gradually decrease during Termination I, F. profunda-based NPP ( Figure 10B) is fairly consistent from LGM to mid-BA, when it starts to decline. The general declining trend is also observed in reconstructed NPP at low resolution at ODP 677B  (Kienast et al., 2006); (E) 230 Th-normalized opal flux (g m -2 yr -1 ) (Kienast et al., 2007); (F) Calculated transfer efficiency (T eff ) following Henson et al. (2012); (G) calculated export production (e f ) using equation 5 (mg C m -2 day -1 ) following Laws et al. (2011) and with 13-point running average smoothing; (H) F. profunda-derived NPP (mg C m -2 day -1 , this work with 13-point running average smoothing). Black dashed lines indicate the late LGM, Termination I and the Holocene. Shaded areas indicate the Heinrich 1 (H1), Older Dryas, Bølling-Allerød (BA) and Younger Dryas (YD). The orange dotted interval highlights high NPP and export production, and low oxygenation of the deep-ocean.
( Figures S2A, C). The NPP-F.profunda in MD02-2529 (8.208°N , 84.125°W, 1619 m depth; Ivanova et al., 2012), however, shows an increase from the LGM to the YD, which then only slowly declines through the Holocene ( Figure S2D). The NPP-F.profunda data, CAR and alkenone concentration at ODP Site 1240 show a significant up swing through the end of H1 and persisting the trough the BA (alkenone accumulation) or YD (NPP-F.profunda and CAR) (Figures 10A-C). Using CAR as a coccolithophore productivity proxy assumes that there is no preservation (carbonate dissolution) and/or dilution bias, in our case supported by high CEX' values ( Figure 8). Both, Anderson et al. (2019) and Jacobel et al. (2020), using different sensitive redox proxies, determined that higher alkenone concentration record prior to the Holocene and deglaciation are the result of changes in bottom water oxygenation at Site 1240. These and other studies in the Pacific Ocean observed enhanced alkenone preservation during poor oxygenation conditions, indicating that alkenone abundance at ODP Site 1240 could be controlled as much by preservation as by primary productivity (Anderson et al., 2019). Although this may account for some of the decline and its dynamics between the LGM and Holocene, the timing of the increase in alkenone accumulation, coincident with other proxies at the end of H1, implies a strong primary signal from increased coccolithophore production at H1. This pattern of coccolithophore productivity increase through H1 to YD is also more consistent with the NPP-F.profunda records from MD02-2529, which peak during the same interval ( Figure S2).
The relatively high primary phytoplankton production from H1 to YD (with estimated NPP values ranging from ca. 1000 to 700 mg C m -2 day -1 , Figure 1H) is consistent with a broad peak in 230 Th normalized organic carbon burial flux at ME0005A-24JC (a site survey core for ODP Site 1240) ( Figure 11D). Even though Kienast et al. (2006) proposed that high palaeoproductivity was fueled by coastal upwelling during those intervals, the records presented here suggest that increased primary production in the EEP phytoplankton community and associated increased organic carbon export was also a strong feature of this interval.
Ba fluxes, an inorganic productivity proxy for carbon rain, have been also measured at ODP Site 1240 (Jacobel et al., 2020) ( Figure 11B b.1), but at a much lower resolution than our NPP record, which makes comparison complicated. At Site MV1014-02-17JC (0.18°S, 85.86°W), however, Ba fluxes show a stepchange at H1 (Figure 11B b.2), more than halving from high LGM values, which broadly correlates with the increase in coccolithophore productivity indices at ODP Site 1240 (NPP-F.profunda, CAR and alkenone accumulation, Figure 10) and a marked rise in Th-normalized opal accumulation, which is another proxy for carbon flux at ME005-24JC ( Figure 11E). Together, these clearly show a switch in productivity and export production regimes, with an increase in calcareous phytoplankton productivity, as well as in opal export, but a marked decline in excess Ba flux. These combined signals are potentially indicative of increased ventilation of sub-thermocline water together with increase nutrient upwelling to the mixedlayer, stimulating phytoplankton production but markedly reducing Ba export. It is worth noting that this is also the start of the decoupling between CAR and NPP-F. profunda (carbonate-based proxies) and alkenone accumulation rates (organic proxy) ( Figure 10) and the start of a decline in 230 Th normalized organic carbon burial flux ( Figure 11D). Again, this all suggests a decoupling between the export efficiency of organic carbon and carbonate through the Termination, whilst primary production remained high. Foraminifera bound d 15 N ( 15 N/ 14 N ratio, a proxy for watercolumn denitrification) measured on Neogloboquadrina dutertrei and N. incompta (average values in Figure 10D) at ME0005A-24JC, show higher values up to the BA suddenly dropping during the YD (Studer et al., 2021). This pattern shows similarities with our CAR and NPP-F.profunda reconstructions, which drop at the end of Termination I (Figures 10A, B). Studer et al. (2021) suggests that planktonic foraminifera d 15 N could be driven by changes in the d 15 N in intermediate water, nitrogen advection and nitrogen utilization by the phytoplankton, on top of the previously inferred d 15 N reduction in water column denitrification during the glacial suggested by Dubois and Kienast (2011) -based on bulk measurements-. The nitrogen utilization by phytoplankton as an explanation for part of the d 15 N variability is supported by our F. profunda-based data, with high steady NPP during Termination I (up to the BA) and a lower values during the Holocene.

Holocene
Almost all of the productivity indicators we considered here are consistent in reconstructing substantially lower Holocene primary productivity and export than that through either the LGM or the Termination, with a long-term trend of decline towards 1.9 ka. There is, however, a broad distinction between those related to organic carbon export and burial (alkenone accumulation, 230 Th normalized organic carbon burial flux, Ba excess flux and Thnormalized opal accumulation) which all decline markedly either through or sharply at the end of Termination I ( Figure 10C and Figures 11B, D, E), and primary productivity proxies based on coccolithophore assemblages or accumulation rates that decrease more steadily through the Holocene (Figures 10A, B). In some instances, in proxies such as our CAR record at ODP Site 1240, there are actual transient increases in the early Holocene, with moderate values persisting through the mid Holocene. By the last 5 kyr, however, even CAR has reached a minimum for the entire record. As with the discussion above for the LGM to Termination interval, this indicates either a decoupling between primary production and organic carbon export efficiency, with the former declining more slowly in the early Holocene, or a potential shift in nutrient regimes that impacted siliceous primary production and organic carbon export at the end of the Termination, but that had less of an impact on the calcareous phytoplankton, which, instead, experienced a more gradual decline in productivity through the Holocene.

Implications for the Export Productivity at ODP Site 1240
In order to better understand the fate of the plankton primary production sinking throughout the water column, we calculated the ratio of export production to total NPP (e f ) at ODP Site 1240 (Equation 5 in this paper) following Laws et al. (2011). For this, we combined our estimates of NPP and SSTa. The reconstructed e f ( Figure 11G) shows a pattern with higher export during the LGM and throughout Termination I, which decreases rapidly at the start of the Holocene, followed by a small increase between 9-5 kyr, to decrease again towards the most recent times. The high glacial e f is followed by a rapid decrease from the BA to the beginning of the Holocene. This rapid decrease in e f is also observed in the transfer efficiency (T eff ) ( Figure 11F), calculated using Th-normalized opal and CaCO 3 fluxes at ME0005-24JC.
We interpret the higher e f values during LGM and the deglaciation ( Figure 11G) as a rise in nutrient supply to the thermocline, as noted by the higher d 15 N in planktonic foraminifera (N. dutertrei and N. incompta) (Studer et al., 2021) (Figure 10D), which resulted in higher NPP and export production (Costa et al., 2017;Jacobel et al., 2020). The lower SST during this interval prevented high remineralization of the organic detritus in the photic zone, leading to a more efficient downward transport of organic detritus and burial in the seafloor (Henson et al., 2012). This higher e f during the LGM and Termination I indicates that the productivity at the EEP contributed actively to decrease the atmospheric CO 2 levels ( Figure 9A), as it was suggested by other studies (e.g., de la Fuente et al., 2017).
The coupled decrease in productivity, e f and T eff during the end of Termination I and the beginning of the Holocene ( Figures 11F, G, H) contrasts with the elevated biogenic opal and organic Carbon fluxes ( Figures 11D, E) at that time. The increased opal productivity between 15-11 kyr has been explained by a higher influence of Si-rich waters from the Southern Ocean transported to the Equator, which stimulated diatom productivity with respect to coccolithophore productivity during the deglaciation (Calvo et al., 2011). It is important to note that the change in ecosystem structure, from calcareous to siliceous plankton, and expected reduction in the carbonate pump (lower calcareous productivity), was not translated into a more efficient export production, as seen by the decrease in T eff ( Figure 11F). A potential explanation for this is the 'packaging effect' (Francois et al., 2002) associated with the shift in the floral community, in which fast-sinking fecal pellets are produced in carbonate-dominated environments, while opal-dominated regions are characterized by very loose aggregates that may sink slower and can be more readily remineralized, reducing the T eff (Henson et al., 2011). Another alternative explanation for the higher opal content, not matching the NPP and e f reconstructions, would be that between 16-9 ka sedimentation rates are higher at site ME0005-24. Faster burial could lead to better preservation of biogenic opal in the sediment (Kienast et al., 2007). Calcium carbonate dissolution has been observed in different Late Pleistocene records of the EEP (e.g., Ivanova et al., 2012;Diz et al., 2018), and would be a mechanism to explain the increase in opal at ME0005A-24JC and decrease in carbonate-based productivity indicators (NPP and CaCO 3 %) between 15-11 kyr. However, good to moderate coccolith preservation has been observed at ODP Site 1240 for the last seven glacial cycles (Loṕez-Otaĺvaro et al., 2008), which is supported by the low CEX' values calculated for our study ( Figure 8J). This suggests that the lysocline was located below 2900 m at Site 1240 for the last ca. 20 kyr. Consequently, we hypothesize that the decoupling between estimates of NPP and e f based on calcareous proxies and % opal is a combined result of shifts in plankton community forced by hydrographic factors, and the better preservation of opal during intervals of higher sedimentation rates in the EEP.
Our records showing higher NPP, e f and T eff during the last part of the LGM and H1 compared to the Holocene, interpreted as a decrease in export production, do not agree with the reconstruction of organic carbon rain at Site 1240 using Thnormalized fluxes of excess of biogenic barium (Ba xs ) by Jacobel et al. (2020), which does not show significant change before and after the LGM ( Figure 11B). To evaluate the signal production and preservation as well as support their interpretation, Jacobel et al. (2020) also measured the authigenic Uranium (aU), a redox proxy, at ODP Site 1240. These authors argue that the lack of correlation between aU and Ba xs ( Figures 11B, C) indicates that the redox proxies are primarily influenced by changes in oxygenation rather than carbon export. However, the Ba xs signal from Site 1240 differs from the higher resolution record at MV1014-02-17JC (Loveley et al., 2017) (Figure 11B), which is around 50 km away and at a similar depth (2846 m) than Site 1240, and shows a more abrupt transition from high to low Ba xs at 15 ka, coinciding with a similar shift from higher to low productivity and export production based on our proxies. We are uncertain about the causes of this difference in Ba xs records from cores located so close, in particular when aU from both sites show the same trends and absolute values ( Figure 11C). It has been argued that under certain oxygen conditions at the sea-floor Ba xs can re-precipitate (McManus et al., 1998), but both sites are have similar deep-ocean oxygen levels, as seen from the similar trend and absolute values in aU. Age model discrepancies are also excluded given the good correlation and absence of lags in both aU records. Although we are not able to explain why the Ba xs records from these sites are different, we can conclude that, given the collective evidence of proxy data at these sites that there was a significant change in surface productivity between the LGM and the Holocene, which also resulted in a change in the flux of carbon to the sea-floor. The reconstruction of bottom ½ CO 2− 3 at Site 1240 (de la Fuente et al., 2017) shows lower values during the LGM and H1 and a rapid increase after BA ( Figure 11A). Accurate comparison is challenging due to the low number of data points in the ½CO 2− 3 reconstruction, but in general lower CaCO 3 and CARs during H1 would be explained by the higher NPP and e f that demand higher levels of carbon respiration at the seafloor, decreasing the ½CO 2− 3 . Therefore, combined productivity and deep carbonate proxies would indicate a more efficient biological pump and enhanced storage of respired carbon in the EEP prior to the deglaciation. Although the YD-Holocene transition is represented by just two points in the ½CO 2− 3 reconstruction, the general trend shows an increase in ½CO 2− 3 , which could be explained by a decrease in productivity as displayed by the NPP reconstruction.

CONCLUSIONS
We have developed a new high-resolution reconstruction of annual sea-surface temperature (SSTa) and net primary productivity (NPP) for the Eastern Equatorial Pacific (EEP) based on novel coccolithophore-based models. The addition of coccolithophore assemblage data from 32 surface sediment samples located in the EEP improves previous SSTcalibrations, and results in higher confidence reconstruction of temperatures warmer than 20°C. Moreover, we observe that the relative abundance of the deep-photic species Florisphaera profunda in these samples fits well with existing calibrations for tropical regions, providing evidence of the robust relationship between this species and satellite-derived NPP within the EEP.
The reconstructed SSTa at ODP Site 1240 using the coccolithbased calibration shows alternating episodes of warming and cooling throughout Termination I, a long-term warming trend across the whole deglaciation and a general slight cooling trend throughout the Holocene. Our SSTa reconstruction differs in the timing and magnitude from other available SST-reconstructions based on biogeochemical and faunal proxies (e.g., foraminiferal Mg/ Ca, foraminiferal assemblages, alkenones). We suggest that these discrepancies could be the result of the proxies differing responses to forcings, seasonal bias and/or preservation artifacts. The reconstructed NPP displays a steady trend during last Glacial and Termination I (up to the Bølling-Allerød) generally decreasing trend to recent times, but differs from other coccolith-based productivity indicators, such as the coccolith accumulation rate (CAR). The increase in coccolithophore productivity proxies at ODP Site 1240 (NPP-F. profunda, CAR and alkenone accumulation) during Heinrich 1, coincides with a marked rise in Th-normalized opal accumulation at ME005-24JC. We attribute this change to increased nutrient upwelling to the mixed-layer caused by the intensification of the trade winds. The calculated ratio of export production to total NPP (e f ) at ODP Site 1240 also indicates higher export during the late Last Glacial Maximum (LGM) and throughout Termination I. The lower SST during this interval prevented high remineralization of the organic detritus in the photic zone, leading to a more efficient downward transport of organic detritus and burial within the sea-floor.
The coupled decrease in reconstructed NPP, e f and transfer efficiency (T eff ) during the end of Termination I, as well as the beginning of the Holocene, contrasts with higher % biogenic opal and organic Carbon fluxes at that time. These changes support the idea of a change in the plankton ecosystem structure at the end of Termination I, from calcareous to siliceous, and a concomitant reduction in the carbonate pump (i.e., lower calcareous productivity). The decrease in export productivity indicators (i.e., Ba xs ) fluxes at the neighboring Site MV1014-02-17JC suggest a decoupling between the surface productivity and its accumulation at the sea-floor, which could be due to a better ventilation of the deep-ocean. This is supported by the redox proxies at ODP Site 1240. All productivity indicators are consistent in reconstructing a lower Holocene primary and exported productivity than during the LGM or the Termination I, consistent with a weakening of wind-driven equatorial upwelling towards recent times.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
The study was designed by J-AF, FJS, MS-P, IH-A, and K-HB. EC analyzed the sediment samples, calculated coccolith abundances and wrote an earlier version of the manuscript. IH-A made the statistical analyses. MS-P and IH-A wrote the manuscript, with input from K-HB, TDJ, J-AF, and FJS. All authors approved the submitted version.