Biological pump variations in the northeastern Indian Ocean since the last glaciation: evidence from coccolithophores

The biological pump (BP) constitutes a key component of the global carbon cycle by mediating the export of organic carbon from the surface ocean to the deep sea. However, long-term variability of BP strength in the northeastern Indian Ocean (NEIO), particularly within the strongly stratified Bay of Bengal (BoB), remains inadequately constrained due to the complex interplay between monsoon forcing and freshwater-induced stratification. Here, we present a high-resolution reconstruction of BP variability over the last ∼ 53 kyr based on primary productivity (PP) estimates derived from quantitative analyses of coccolithophore assemblages in a sediment core from the southern BoB. Coccoliths are abundant and well preserved throughout the record, with assemblages dominated by Florisphaera profunda, Gephyrocapsa oceanica, and Emiliania huxleyi. The reconstructed PP record exhibits pronounced glacial–interglacial (G–IG) variability, characterized by reduced productivity during the Last Glacial Maximum, followed by a progressive increase during the deglaciation and a peak in the Early Holocene. Principal component analysis of productivity records across the Indian Ocean (IO) identifies G–IG climate variability as the dominant mode of basin-scale productivity change while distinguishing secondary atmospheric and hydrological controls on regional productivity. Comparison with PP records from other sectors of the IO further highlights marked spatial heterogeneity in BP responses to basin-scale climatic forcing. In the NEIO, productivity was primarily controlled by G–IG changes in upper-ocean stratification and nutrient availability driven by sea level fluctuations and freshwater forcing, whereas secondary variability reflects the distinct atmospheric and hydrological expressions of the Indian Summer Monsoon. Our findings indicate that regional hydrographic processes governed the expression of basin-scale climatic forcing, distinguishing the freshwater-controlled NEIO from the wind-driven upwelling systems of the western tropical IO.

1 Introduction

The marine biological pump (BP) constitutes a fundamental component of the Earth's carbon cycle, governing the exchange of carbon between the ocean and atmosphere through the vertical transport of organic matter from surface waters to the deep ocean and underlying sediments. The efficiency of this process is closely coupled to marine primary productivity (PP), which determines the production and subsequent export of organic carbon to depth. As the primary drivers of marine photosynthesis, phytoplankton exert a fundamental control on global carbon cycling by modulating BP intensity (Litchman et al., 2015). Variations in BP strength can therefore alter marine ecosystem structure and significantly influence the ocean's capacity for atmospheric CO₂ sequestration, whereas diminished BP activity limits long-term carbon storage and weakens climate regulation (Häder et al., 2014). Despite its importance, the long-term evolution of BP strength and its governing mechanisms remain insufficiently constrained in several climatically sensitive regions. In particular, it remains unclear how large-scale climatic forcing and regional oceanographic processes interact to regulate BP variability over glacial–interglacial (G–IG) timescales. Addressing this gap is essential for resolving the links between ocean productivity, carbon sequestration, and climate change (Séférian et al., 2014).

On G–IG timescales, BP variability reflects the combined influence of physical, chemical, and biological processes that regulate upper-ocean structure and nutrient availability. Climate-driven changes in sea level, wind stress, and hydrological balance reorganize ocean circulation, stratification, and vertical nutrient transport, thereby modulating PP and associated carbon export (Grant et al., 2014; Farmer et al., 2021; Letscher et al., 2023; Shankle et al., 2025). In tropical regions, monsoon systems exert a particularly strong control by altering surface stratification and nutrient supply through coupled variations in wind forcing and freshwater input (Le Mézo et al., 2017; Windler et al., 2021). The interaction of these drivers produces spatially heterogeneous BP responses, especially between upwelling-dominated margins and strongly stratified basins (Crosta and Shemesh, 2002; Le Mézo et al., 2017; Letscher et al., 2023;). However, the relative importance of these mechanisms during past climate transitions remains incompletely resolved, limiting our ability to disentangle the processes linking ocean productivity, atmospheric CO₂ variability, and last glaciation (LG) climate evolution.

The northeastern Indian Ocean (NEIO) constitutes a climatically sensitive setting where BP variability is strongly modulated by monsoon-driven atmospheric circulation and pronounced salinity stratification arising from substantial freshwater input, particularly within the Bay of Bengal (BoB) (Shetye and Shenoi, 1988; Kumar et al., 1996; Priya et al., 2015; Gordon et al., 2016). In contrast to upwelling-dominated sectors of the Indian Ocean (IO), where wind-driven Ekman transport enhances nutrient supply and sustains elevated biological production, the NEIO is characterized by a shallow mixed layer and a deep nutricline that generally suppress nutrient entrainment into the euphotic zone, thereby limiting BP strength (Prasanna Kumar et al., 2002; Wiggert et al., 2005; Phillips et al., 2014; Da Silva et al., 2017; Bolton et al., 2022). Paleoceanographic investigations from this region document diverse G–IG productivity patterns; however, the extent to which BP variability in the BoB reflects basin-scale climatic forcing versus regional hydrographic controls remains poorly constrained (Madhupratap et al., 2003; Prasad, 2004; Prasanna Kumar et al., 2004; Gauns et al., 2005; Sarma et al., 2019; Zhou et al., 2020). This ambiguity underscores the need for high-resolution, proxy-based reconstructions capable of resolving the mechanisms governing BP variability in this strongly stratified, monsoon-influenced basin.

Coccolithophores, as calcifying phytoplankton, play a dual role in marine carbon cycling by contributing to both organic carbon fixation and calcium carbonate production (Balch et al., 2018; Rost and Riebesell, 2004). Through their involvement in the organic and carbonate pumps, they facilitate the transfer of carbon from surface waters to the deep ocean (Srivastava et al., 2025). Their calcareous plates (coccoliths) are well preserved in marine sediments, providing robust proxies for reconstructing past oceanographic and biogeochemical conditions (Young, 1998; Jacques and Luc, 2007). Variations in coccolithophore assemblage composition and fluxes are closely linked to changes in upper-ocean stratification and nutrient dynamics, offering valuable insights into BP variability. In particular, Florisphaera profunda serves as a key indicator species, with its relative abundance reflecting nutricline depth and providing a quantitative basis for reconstructing PP and associated BP intensity (Molfino and McIntyre, 1990; Beaufort et al., 1997; Stoll et al., 2007; Zhang et al., 2016; Hernández-Almeida et al., 2019).

This study reconstructs BP variability in the NEIO since the LG (∼ 53 kyr) using coccolithophore assemblages and quantitative PP estimates derived from the high-resolution sediment core YDY05. Integration of this record with productivity reconstructions from other sectors of the IO enables evaluation of the relative roles of basin-scale climatic forcing and regional hydrographic processes in regulating BP dynamics within the BoB. Furthermore, this approach facilitates identification of the mechanisms underlying the contrasting productivity responses between the southern BoB and upwelling-dominated regions of the IO.

2 Modern oceanic condition

The BoB, located in the NEIO, is a U-shaped, semi-enclosed basin bounded by land on three sides and connected to the open IO along its southern margin (Fig. 1). The modern BP in this region is primarily governed by monsoon-driven circulation and its interaction with freshwater forcing, upper-ocean stratification, and nutrient redistribution, which collectively regulate primary production and particulate organic carbon (POC) export (Thushara et al., 2019).

Figure 1 (a) Geographic setting and bathymetric map of the Indian monsoon climate zone, including the Bay of Bengal (BoB), Andaman Sea, and Arabian Sea (AS) of the Indian Ocean (IO). The map is created by Arc GIS software. The location of the sediment core YDY05 is marked by a red star, and additional comparable published cores (SK157-16, MGS22/01, YDY09, BoB-24, MD77-176, MD161-19, SO188-17286-1, MD90-0963, MD85-668, GeoB12613-1, and SO139-74KL) are highlighted by a brown circle. White and black arrows show surface ocean circulation during the summer monsoon and winter monsoon. The red arrows show the eastward flow of waters from the AS into the BoB via the Somali Current (SC), West Indian Coastal Current (WICC), and Southwest Monsoon Current (SMC). The dashed black line indicates the direction of sea surface currents: the East Indian Coastal Current (EICC), East African Coastal Current (EACC), Northeast Madagascar Current (NEMC), South Equatorial Current (SEC), South Equatorial Countercurrent (SECC), and Indonesian Throughflow (ITF). (b–i) The modern climatology of the IO. (b, c) Mean precipitation rate for Northern Hemisphere (NH) summer (July–August–September) and winter (January–February–March), respectively. Data are from NCEP-DOE Reanalysis 2 (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis2.html, last access: 8 February 2026). (d, e) Mean surface wind speed rate for NH summer and winter, respectively. Data are from the CPC Merged Analysis of Precipitation (https://psl.noaa.gov/data/gridded/data.cmap.html, last access: 8 February 2026). Average monthly phosphate (f, g) and sea surface salinity (h, i) of summer and winter, respectively. These climatological datasets are compiled from the World Ocean Atlas 2023 and visualized using Ocean Data View 4 © (Schlitzer, 2016; https://odv.awi.de, last access: 29 June 2025).

The seasonal monsoon exerts a fundamental control on circulation and nutrient distribution in the BoB. During the Indian Summer Monsoon (ISM), strong southwesterly (SW) winds reorganize surface circulation (Schott et al., 2009), while phosphate concentrations increase markedly in the western BoB (Fig. 1), indicating enhanced nutrient availability associated with wind-driven upwelling and mixing. Consequently, this region supports relatively high PP despite the presence of surface stratification. In contrast, phosphate concentrations at Site YDY05 are lower than those in the western BoB but remain higher than in much of the central basin, suggesting that nutrient supply at the study site is maintained primarily through monsoon-driven circulation rather than local upwelling. During the winter monsoon, northeasterly (NE) winds promote upper-ocean mixing, providing an additional seasonal mechanism for nutrient replenishment and sustaining PP and POC export (Shetye et al., 1996).

In addition to wind-driven circulation, salinity is a key regulator of nutrient availability in the BoB. During the summer monsoon, intense precipitation (∼ 20 mm d⁻¹), together with discharge from the Ganges–Brahmaputra–Meghna and Irrawaddy river systems, substantially lowers sea surface salinity compared with the Arabian Sea (AS) (∼ 1–6 mm d⁻¹; Fig. 1b–c; Randel and Park, 2006; Sarma et al., 2016; Thushara and Vinayachandran, 2016). The resulting low-density surface layer overlies relatively saline subsurface waters, forming a strong halocline layer and barrier layer that inhibit vertical mixing. Consequently, nutrient-rich subsurface waters remain isolated from the euphotic zone, limiting phytoplankton growth and reducing PP despite favorable light conditions. When freshwater input decreases, surface salinity increases, upper-ocean stratification weakens, and nutrient entrainment into the euphotic zone becomes more efficient, promoting higher productivity. Modern observations illustrate this relationship: during the spring intermonsoon, PP ranges from approximately 13.2 to 173.8 mmol C m⁻² d⁻¹, whereas POC export remains relatively low (0–7.7 mmol C m⁻² d⁻¹), reflecting persistent nutrient limitation associated with upper-ocean stratification (Subha Anand et al., 2017). Consequently, the BoB generally exhibits lower productivity than the AS despite experiencing similar monsoon forcing (Prasanna Kumar et al., 2002).

Although freshwater-induced stratification suppresses vertical nutrient supply across much of the basin, the southern BoB is less directly influenced by river discharge than the northern BoB. Instead, nutrient availability at Site YDY05 is strongly enhanced by lateral advection and monsoon-driven circulation. During the ISM, the Southwest Monsoon Current (SMC), an eastward extension of the Somali Current, flows south of Sri Lanka and transports nutrient-rich AS waters into the southern BoB at velocities exceeding 40 cm s⁻¹ and transport rates of approximately 8–15 Sv (Schott et al., 1994; Vinayachandran et al., 2004). Mesoscale eddies associated with the SMC further inject nutrients into the upper ocean (Prasanna Kumar et al., 2004; McCreary et al., 2009), while exchange with the equatorial IO provides an additional nutrient source (Sanchez-Franks et al., 2019). These processes help offset the suppressing effect of freshwater stratification and explain the relatively elevated nutrient availability observed at Site YDY05 during the summer monsoon (Fig. 1). Thus, BP variability in the southern BoB reflects the balance between monsoon-driven nutrient supply and freshwater-controlled salinity stratification.

3 Materials and methods

3.1 Sediment core and age model

This study is based on quantitative analyses of coccolithophore assemblages from piston core YDY05, retrieved from the southern BoB (9.99° N, 90.32° E; water depth: 3300 m; Fig. 1a) during an open research cruise in the eastern IO conducted by R/V Shiyan 3, affiliated with the South China Sea Institute of Oceanology, Chinese Academy of Sciences, China. The 2.39 m long core is composed predominantly of olive-gray clay and silty clay with relatively high carbonate content. The core was continuously subsampled at 2 cm intervals, yielding a total of 120 samples.

The chronological framework of core YDY05 is constrained by six accelerator mass spectrometry (AMS) radiocarbon dates obtained from mixed planktonic foraminifera (Table S1 in the Supplement; Luo et al., 2018; Devendra et al., 2019). Radiocarbon ages were calibrated to calendar years using the Marine20 calibration dataset implemented in CALIB version 8.20, applying a global marine reservoir correction of 51 years and a regional ΔR of −61 years to account for spatial variability in surface ocean radiocarbon concentrations (Dutta et al., 2001). An age–depth model was constructed using the Bayesian software Bacon within the R statistical environment (version 4.3.1; Fig. 2), which integrates calibrated ages with prior constraints on sediment accumulation variability (Blaauw and Christen, 2011). Marine Isotope Stage (MIS) boundaries in this study are defined following the LR04 benthic δ¹⁸O stack of Lisiecki and Raymo (2005). Accordingly, MIS 1 is defined as 0–14 kyr, MIS 2 as 14–29 kyr, and MIS 3 as 29–57 kyr. These boundaries are applied consistently across all figures and age models in this study.

Figure 2 The age–depth model. The model and sedimentation rates of core YDY05 are established based on the data obtained from Bayesian analysis with the Bacon software. The dotted-dash red line represents the mean, and the solid gray line indicates the median age. Besides this, the model's 95 % probability interval is specified by the area between black and green lines, showing the maximum and minimum ages.

3.2 Coccolith analysis and PP reconstructions

Quantitative coccolith analyses were performed following the standard drop technique for slide preparation (Bordiga et al., 2015). Approximately 15 mg of oven-dried sediment was weighed and transferred into a 15 mL centrifuge tube. To disaggregate particles and remove organic matter, 10 mL of a 2 % ammonia solution (pH ≈ 10) and 300 µL of 10 % hydrogen peroxide were added. The suspension was ultrasonicated for 15 s and subsequently allowed to settle for 24 h. A 300 µL aliquot from the upper suspension was pipetted onto a 24 × 24 mm² coverslip, evenly distributed, and dried at 40 °C. The coverslip was then permanently mounted onto a glass slide using Norland Optical Adhesive (NOA61).

Coccolith identification and counting were conducted under an Olympus BX53 polarizing-light microscope at 1000× magnification using both cross-polarized and transmitted light. Taxonomic identification followed Young (1998) and Young et al. (2003) as compiled in the Nannotax3 online database (Nannotax3 website, 2025). A minimum of 300 coccoliths were counted per sample to ensure statistically robust representation of assemblage composition. Absolute coccolith abundance (specimens g⁻¹ dry sediment) and relative abundances (%) of individual taxa were calculated for each sample to evaluate temporal variations in assemblage structure.

Absolute coccolith abundance (AA) was calculated using the following equation:

$$\begin{array}{}\text{(1)}& \mathrm{AA}={\displaystyle \frac{N\times S\times V}{s\times v\times m}},\end{array}$$

where N is the number of coccoliths counted, S is the total coverslip area (24 × 24 mm²), V is the total suspension volume (10.3 mL), s is the cumulative area of the counted fields of view, v is the volume of suspension dispensed onto the coverslip (300 µL), and m is the dry sediment mass (g). Although all identifiable coccolith taxa were recorded, subsequent analyses focus on the three ecologically significant and numerically dominant species F. profunda, Emiliania huxleyi, and Gephyrocapsa oceanica. Relative abundance (RA) was calculated as follows:

$$\begin{array}{}\text{(2)}& \mathrm{RA}={\displaystyle \frac{ni}{N}}\times \mathrm{100},\end{array}$$

where n_i is the number of coccoliths of species i, and N is the total number of coccoliths counted in the sample. Coccolith fluxes (C_f) were calculated to quantify temporal variations in coccolithophore export production. Fluxes (coccoliths cm⁻² kyr⁻¹) were derived by multiplying total absolute coccolith abundance (A) by the mass accumulation rate (MAR):

$$\begin{array}{}\text{(3)}& C_{\mathrm{f}}=A\times \mathrm{MAR}.\end{array}$$

Mass accumulation rate was calculated as follows:

$$\begin{array}{}\text{(4)}& \mathrm{MAR}=\mathrm{LSR}\times \mathrm{DBD},\end{array}$$

where LSR is the linear sedimentation rate obtained from the age–depth model, and DBD is dry bulk density. As direct dry bulk density measurements were unavailable for core YDY05, an average value of 1.27 g cm⁻³ derived from a nearby BoB core with comparable lithology and depositional conditions was applied (Khan et al., 2025). While this assumption may influence absolute flux magnitudes, it does not affect relative downcore trends, which constitute the primary focus of this study.

PP was reconstructed from the relative abundance of F. profunda (Fp %), a proxy indicative of nutricline depth and surface productivity in tropical oceans. The empirical relationship between Fp % and PP was first established for the AS (Beaufort et al., 1997). In this study, PP over the last ∼ 53 kyr was estimated using an updated Fp %–PP calibration for the tropical IO, derived from a global core top dataset and constrained by satellite-based productivity estimates (Hernández-Almeida et al., 2019). PP was calculated using the following equation:

$$\begin{array}{}\text{(5)}& \mathrm{PP}=\left[{\mathrm{10}}^{\left(\mathrm{3.27}-\mathrm{0.01}\times \mathrm{Fp}\mathit{\%}\right)}\right]\times {\displaystyle \frac{\mathrm{365}}{\mathrm{1000}}},\end{array}$$

where PP is expressed in g C m⁻² yr⁻¹. The calibration exhibits a strong correlation between estimated and observed productivity (r²=0.74) with a residual standard error of 1.5 g C m⁻² yr⁻¹ (Hernández-Almeida et al., 2019).

3.3 Statistical analysis

Principal component analysis (PCA) was applied to Late Quaternary PP proxy records from six sediment cores distributed across the IO, including YDY05, MGS22/01, SK157/16, MD161-19, GeoB12613-1, and SO139-74KL (Table S2). The productivity proxies comprised coccolith-driven PP, biogenic barium (Ba_bio), total organic carbon (TOC) accumulation rate, and coccolith fraction (CF) $\mathrm{Sr}/\mathrm{Ca}$. Published records from MD77-176, MD90-0963, and MD85-668 were excluded because MD77-176 covers only the last ∼ 26 kyr, numerical data for MD90-0963 were unavailable, and MD85-668 has a substantially lower temporal resolution than the other records.

To minimize the influence of different temporal resolutions among the proxy records, all datasets were linearly interpolated to a common temporal interval of 1.7 kyr (1700 years) using OriginPro (OriginLab Corporation, Northampton, MA, USA) prior to PCA. The interpolated datasets were subsequently standardized by subtracting the mean and dividing by the standard deviation to remove scale-dependent effects and ensure equal weighting among variables (Sokal and Rohlf, 1995). PCA was then performed using PAST version 4.03 (Hammer et al., 2001). Principal component scores and loadings were used to identify the dominant modes of productivity variability and their relationships among the different proxy records. To facilitate interpretation of the principal components, a second PCA was performed using the PC1 and PC2 scores together with selected climatic and oceanographic proxy records. The resulting PCA biplot was used to evaluate the relationships between the principal components and the major climatic and oceanographic drivers.

4 Results

4.1 Chronology and sedimentation rates

Linear extrapolation suggests a basal age of ∼ 52.68 kyr, indicating that the core preserves a continuous late Quaternary depositional record. Sedimentation rates vary throughout the sequence, ranging from ∼ 3.7 to ∼ 8.6 cm kyr⁻¹. Elevated rates (∼ 8.6–5.1 cm kyr⁻¹) occur between ∼ 27 and 16 kyr, whereas lower rates (∼ 3.7 cm kyr⁻¹) characterize the Holocene interval (∼ 11–5 kyr) (Fig. 2).

4.2 Coccolith absolute abundances and flux

Coccoliths are abundant and well preserved in all analyzed samples from core YDY05. Total coccolith concentrations range from 3.72 × 10⁹ to 16.7 × 10⁹ coccoliths g⁻¹, with a mean value of 9.55 × 10⁹ coccoliths g⁻¹ (Fig. 3a). The assemblages are dominated by F. profunda, E. huxleyi, and G. oceanica, which, together, account for ∼ 85 %–97 % of the total coccolith community. Among these taxa, F. profunda exhibits the highest mean concentration (5.56 × 10⁹ coccoliths g⁻¹), exceeding those of E. huxleyi and G. oceanica (Fig. S1 in the Supplement). Species-specific trends indicate that F. profunda predominates prior to the Last Glacial Maximum (LGM), followed by a marked decline to 1.26 × 10⁹ coccoliths g⁻¹ thereafter. In contrast, E. huxleyi and G. oceanica show progressive increases, with pronounced maxima at ∼ 11 kyr. Total coccolith fluxes range from 9.5 × 10⁹ to 4.2 × 10¹⁰ coccoliths cm⁻² kyr⁻¹, with an average of 2 × 10¹⁰ coccoliths cm⁻² kyr⁻¹ (Fig. 3d). Fluxes are lowest during ∼ 53–27 kyr and reach peak values between ∼ 18 and ∼ 11 kyr.

Figure 3 (a) Absolute abundance of the total no. of coccolith species. (b) Absolute abundance of Noelaerhabdaceae species (Gephyrocapsa oceanica, Emiliania huxleyi). (c) Relative abundance (%) of Florisphaera profunda. (d) Coccolith flux. (e) Primary productivity (PP) records. The curves (a, b, d) are smoothed using the LOESS (locally estimated scatterplot smoothing) algorithm implemented in PAST software, with a smoothing factor (span) of 0.1. The smoothing factor indicates that each local regression was fitted using approximately 10 % of the neighboring data points. The smoothed curves are shown for visual interpretation only. The color bars show millennial-scale intervals: Deglaciation and Last Glacial Maximum (LGM).

4.3 Reconstructed PP of indicative species

PP was reconstructed from the relative abundance of the dominant species F. profunda. This taxon constitutes, on average, ∼ 57 % of the assemblages, whereas E. huxleyi and G. oceanica generally contribute less than 25 % (Fig. S2 in the Supplement). The most pronounced compositional changes occur between ∼ 34 and ∼ 11 kyr, particularly during ∼ 27–25 kyr, when F. profunda increases sharply to ∼ 80 %, accompanied by a decline in E. huxleyi to ∼ 4 % and a slight decrease in G. oceanica from ∼ 18 % to ∼ 14 %. Based on these variations, three principal intervals are identified: (i) between ∼ 52 and ∼ 24 kyr, F. profunda predominates, with a mean relative abundance of ∼ 63 %. This interval is characterized by generally stable values around ∼ 60 %, interrupted by two prolonged intervals of reduced abundance at ∼ 49–45 kyr and ∼ 35–32 kyr, and a pronounced maximum during ∼ 27–25 kyr. (ii) The interval from ∼ 20 to ∼ 11 kyr is characterized by pronounced variability among all three taxa, with amplitudes of ∼ 22 % for F. profunda, ∼ 25 % for E. huxleyi, and up to ∼ 46 % for G. oceanica. (iii) Between ∼ 11 and ∼ 4 kyr, F. profunda increases progressively, reaching ∼ 45 % at ∼ 8 kyr, followed by a gradual decline toward ∼ 4 kyr.

The reconstructed PP record exhibits clear G–IG variability (Fig. 3). Values range from ∼ 108 to ∼ 412 g C m⁻² yr⁻¹, with minimum productivity during the LGM and peak values during the Early Holocene (EH). Between ∼ 53 and ∼ 27 kyr, PP averages ∼ 166 g C m⁻² yr⁻¹, followed by a decline to ∼ 108 g C m⁻² yr⁻¹ at ∼ 26–25 kyr. A marked increase occurs during the deglaciation, with PP rising from ∼ 146 g C m⁻² yr⁻¹ at ∼ 18 kyr to ∼ 412 g C m⁻² yr⁻¹ at ∼ 11 kyr. Although productivity remains relatively elevated during the EH, it subsequently declines to ∼ 239 g C m⁻² yr⁻¹ during the mid-Holocene.

4.4 PCA

PCA was applied to standardized productivity-related proxy records from six sediment cores to identify the dominant modes of productivity variability across the IO. The first two principal components explain 80 % of the total variance, with PC1 and PC2 accounting for 67 % and 13 %, respectively. PC1 exhibits a pronounced inverse relationship with the global benthic δ¹⁸O stack, whereas PC2 broadly follows the temporal evolution of the Indian Monsoon Stack (Fig. 4A, B). The spatial distribution of PCA loadings reveals distinct regional patterns among the productivity records (Fig. 4C, D). PC1 is characterized by positive loadings for the NEIO records (YDY05, MGS22/01, SK157/16, and MD161-19) and negative loadings for the western and eastern tropical IO records (GeoB12613-1 and SO139-74KL). In contrast, PC2 exhibits the strongest positive loading for GeoB12613-1; moderate positive loadings for MGS22/01, SK157/16, and MD161-19; a weak positive loading for YDY05; and a near-zero loading for SO139-74KL.

Figure 4 PCA results of primary productivity proxy records from the Indian Ocean. (a) PC 1 (dark green) and global benthic stack of core LR04 (orange) (Lisiecki and Raymo, 2005). (b) PC 2 (green) and the Indian summer monsoon stack (orange) (Clemens and Prell, 2003). (c) Loading values of PC1 and (d) PC2 for northeastern Indian Ocean (NEIO), western tropical Indian Ocean (WIO), and eastern tropical Indian Ocean (EIO) sites. The dark-green and orange color sites show positive and negative loadings, respectively.

5 Discussion

5.1 Coccolithophore-derived BP strength variations since 53 kyr

Coccolithophore assemblages from core YDY05 provide a robust basis for reconstructing BP variability in the NEIO. The observed species composition is consistent with previous studies from the IO, where F. profunda, G. oceanica, and E. huxleyi dominate (Beaufort et al., 1997; Rogalla and Andruleit, 2005; Andruleit et al., 2008; Tangunan et al., 2017; Zhou et al., 2020). Unlike commonly used productivity proxies such as TOC, Ba(bio), BSi, carbonate accumulation rates, and planktonic foraminiferal δ¹³C, which may be influenced by preservation, dissolution, diagenesis, or terrigenous inputs (Mackensen and Bickert, 1999; Conley and Schelske, 2002; Galy et al., 2007; Honjo et al., 2014; Liguori et al., 2016), coccolith assemblages are particularly well suited for reconstructing BP variability because changes in species composition directly reflect variations in nutrient availability, nutricline depth, and upper-ocean stratification. This ecological sensitivity is especially valuable in the BoB, where monsoon-driven stratification is the primary regulator of nutrient supply and PP. In particular, the deep-dwelling species F. profunda serves as a key indicator of nutricline depth and stratification, with higher relative abundances reflecting reduced nutrient supply to the euphotic zone (Baumann et al., 2005; Hernández-Almeida et al., 2019). Its persistently high abundance therefore suggests that strongly stratified conditions prevailed in the NEIO over the last ∼ 53 kyr.

Intervals of elevated F. profunda are interpreted to reflect intensified stratification and nutrient limitation, resulting in lower reconstructed PP. In contrast, enhanced productivity corresponds to increased contributions from upper-euphotic Noelaerhabdaceae taxa (i.e., G. oceanica and E. huxleyi) (Fig. 3), which preferentially thrive under more nutrient-enriched surface conditions (Bolton et al., 2010). Moreover, the coherent covariation among total coccolith abundance, flux, and reconstructed PP indicates that assemblage dynamics primarily reflect variations in upper-ocean structure rather than preservation effects or sedimentary dilution. The elevated coccolith abundance during the LGM, largely driven by F. profunda, therefore reflects enhanced upper-ocean stability rather than increased productivity.

The reconstructed PP record exhibits a pronounced G–IG pattern supported by independent proxies. Productivity is relatively elevated during MIS 3 and MIS 1 but declines markedly during MIS 2, reaching minimum levels during the LGM. A progressive increase during the deglaciation culminates in peak values during the EH. Throughout this interval, both the F. profunda and PP records exhibit relatively smooth variations, without clearly expressed oscillations corresponding to the Bølling–Allerød or Heinrich Stadial 1 signals. This pattern is consistent with the BSi record from the same core (Fig. 5b). A short-lived decline in PP between ∼ 26 and 25 kyr (∼ 10 %–15 % lower than the preceding interval) is also observed. The cause of this brief reduction remains uncertain and is therefore not interpreted further.

Figure 5 (a) Reconstructed primary productivity (PP) of this study. (b) Biogenic Silica of YDY05 (Zhang et al., 2022). (c) Total organic carbon (TOC) of core MD161-19 (Da Silva et al., 2017). (d) Biogenic barium (Ba_bio) records of core MGS 22/01 (Fernandes et al., 2025). (e) Biogenic barium (Ba_bio) records of core SK157-16 (Banerjee et al., 2024). (f) PP records of core MD90-0963 (Beaufort et al., 1997). (g) PP records of core MD77-176 (Zhou et al., 2020). (h) PP records of core GeoB12613-1 (Tangunan et al., 2017). (i) PP records of core MD85-668 (Rickaby et al., 2007). (j) PP records of core SO139-74 (Andruleit et al., 2008). The dark-green and red color curves show Glacial–Holocene PP increasing and decreasing trend, respectively. The color bars show millennial-scale intervals: Deglaciation and Last Glacial Maximum (LGM).

The magnitude and pattern of reconstructed PP further support the regional robustness of the proxy. The mean PP (∼ 192 g C m⁻² yr⁻¹) closely matches modern satellite-derived estimates for the southern BoB (Saxena et al., 2020). The EH maximum coincides with increased coccolith abundance and flux, indicating enhanced surface production and export potential, followed by a gradual decline toward the mid-Holocene. Collectively, these results demonstrate a coherent G–IG evolution of BP strength in the NEIO, primarily governed by changes in upper-ocean stratification and nutrient availability.

5.2 BP strength variability across the IO

BP variability across the IO exhibits pronounced spatial heterogeneity, reflecting differences in basin geometry, regional circulation, freshwater flux, and monsoon-driven hydrographic structure (Kumar et al., 1996; Seo et al., 2009; Fournier et al., 2017). To contextualize productivity changes at Site YDY05, we compared our reconstructed PP record with published records from different sectors of the basin (Fig. 5; Table S2). Within the NEIO, the TOC MAR record from core MD161-19, together with the Ba_bio records from core MGS 22/01 and SK157-16, indicates elevated productivity during MIS 1 and MIS 3 and reduced productivity during MIS 2, consistently with the pattern observed in this present study. Nevertheless, spatial deviations occur among BoB records, reflecting contrasting regional controls on nutrient supply. For example, the northern BoB core MD77-176 exhibits relatively stable productivity during MIS 2. According to Zhou et al. (2020), this pattern was influenced by lowered sea level during the LGM, which shifted the mouths of the Irrawaddy and Salween rivers closer to the core site, thereby enhancing river-derived nutrient delivery during episodic increases in discharge. As sea level rose during the deglaciation, the direct effect of riverine nutrients diminished, and productivity became increasingly controlled by freshwater-induced salinity stratification. In contrast, core MD161-19 from the northern BoB records enhanced productivity during glacial periods primarily because reduced monsoonal freshwater input weakened the low-salinity surface layer, promoted vertical mixing and nutrient entrainment, and thereby increased biological productivity (Da Silva et al., 2017). These contrasting records highlight that productivity in the BoB is strongly modulated by regional hydrographic conditions, including riverine nutrient supply, salinity stratification, and shelf configuration, in addition to basin-scale climatic forcing.

Beyond the NEIO, productivity evolution shows contrasting regional expressions. Records from the equatorial IO (MD90-0963) display a G–IG pattern broadly comparable to that observed at Site YDY05, with reduced productivity during the LGM and enhanced productivity during the deglaciation and Holocene, indicating coherent low-latitude variability across the basin. In contrast, upwelling-dominated records from the western tropical and eastern IO (GeoB12613-1, MD85-668, and SO139-74KL) frequently exhibit the opposite behavior, characterized by elevated productivity during glacial intervals and comparatively reduced values during interglacial periods. These divergent patterns highlight pronounced spatial heterogeneity in PP and associated BP strength across the IO and emphasize the influence of regional hydrographic regimes on their basin-scale variability (DiNezio et al., 2018; Vázquez et al., 2024).

Figure 6 Influencing factors for biological pump strength. (a) PC1 curve. (b) Global benthic level (Lisiecki and Raymo, 2005). (c) Global sea level (Miller et al., 2020). (d) PC2 curve. (e) Indian summer monsoon stack (Clemens and Prell, 2003). (f) δ¹⁸Osw-ivc data of core SO188-17286-1 (Lauterbach et al., 2020). (g) δ18Osw data of core BoB-24 (Liu et al., 2021). (h) Records of speleothem δ¹⁸O from Sanbao Cave in central China (Cheng et al., 2016). (i) Primary productivity of YDY05 (this study). Different color bars show the most obvious productivity transformation among these proxy records from different millennial-scale intervals: Mid-Holocene (MH), Early Holocene (EH), Heinrich Stadial 1 (HS-1), Bølling-Allerød (B/A), the Younger Dryas (YD), and Last Glacial Maximum (LGM).

Figure 7 Principal component analysis biplot showing the relationships among productivity- and climate-related proxy records used in this study. Black circles represent sample scores; blue and red arrows represent the loading vectors of the proxy variables. The direction and length of each loading vector indicate the correlation and relative contribution of the corresponding variable to the first two principal components.

5.3 Controls on BP variability since the LG

The long-term PP variations reconstructed from core YDY05 primarily reflect changes in upper-ocean hydrography that regulated nutrient availability in the southern BoB. In this region, the vertical exchange of nutrients is strongly controlled by salinity-driven stratification, which determines mixed-layer depth and the position of the thermocline and nutricline (Thushara et al., 2019). Consequently, even relatively small changes in freshwater balance and water column stability can substantially alter nutrient delivery to the euphotic zone, thereby affecting BP.

Figure 8 Conceptual model illustrating the influence of monsoon variability on BP strength and carbonate sedimentation in the northeastern Indian Ocean. Glacial conditions with weakened summer monsoon circulation correspond to reduced primary productivity and weaker BP strength, whereas stronger summer monsoon conditions during the Holocene promoted enhanced productivity, intensified carbon export, and increased BP strength.

During glacial intervals, the expansion of continental ice sheets and an approximately 120 m decline in global sea level (Fig. 6c; Miller et al., 2020) fundamentally reorganized Indo-Pacific ocean circulation. A reduced sea level restricted the Indonesian Throughflow (ITF), decreasing the transport of relatively fresh Pacific surface waters into the IO and increasing basin-scale surface salinity (Nuber et al., 2023). However, hydrographic conditions in the southern BoB were governed primarily by regional freshwater balance rather than ITF variability alone. Reduced precipitation, river discharge, and shelf–basin connectivity maintained a strongly stratified upper ocean, deepened the nutricline, and limited nutrient exchange with surface waters (Da Silva et al., 2017; Galy et al., 2008; Phillips et al., 2014). Together, these processes suppressed phytoplankton production and organic carbon export, consistently with the weakened BP reconstructed during the LGM.

Deglaciation initiated a progressive reorganization of regional hydrography. Rising sea level restored shelf–basin connectivity, while cyclonic eddies, island-induced upwelling around the Maldives, and equatorial circulation associated with the Wyrtki Jet periodically shoaled the thermocline and enhanced nutrient delivery to surface waters (Chen et al., 2016; Thushara et al., 2019). These processes promoted higher PP and export production during the EH, strengthening the BP through increased sequestration of POC (Honjo et al., 2014). However, productivity gradually declined after the EH despite intensified summer monsoon conditions. Increased monsoonal precipitation enhanced freshwater input and reinforced upper-ocean stratification, offsetting the positive effects of stronger monsoon circulation by suppressing vertical nutrient transport (Da Silva et al., 2017; Thushara et al., 2019). This indicates that BP variability in the southern BoB was ultimately governed by the balance between nutrient supply and freshwater-induced stratification rather than by monsoon intensity alone.

To place the hydrographic controls identified at Site YDY05 within a basin-wide context, we compared productivity records from across the IO using PCA. As the productivity records differ in terms of their original temporal resolution, the PCA is interpreted primarily in terms of coherent basin-scale and orbital-scale variability. Accordingly, although interpolation provides a common temporal framework for comparison, higher-frequency fluctuations should be interpreted with appropriate caution, particularly for records with coarser original sampling resolution (please see details in Sect. 3.3). The dominant pattern of variability (PC1) closely tracks the global benthic δ¹⁸O stack (Fig. 6a, b), indicating that G–IG climate change exerted the primary control on long-term productivity across the basin. However, the contrasting spatial distribution of PC1 loadings reveals that this common climatic forcing produced regionally distinct hydrographic responses. Productivity records from the NEIO exhibit positive loadings, whereas those from the western and eastern tropical IO display negative loadings (Fig. 4B), indicating fundamentally different mechanisms regulating nutrient supply and biological productivity.

Although the western and eastern tropical IO differ in terms of their regional oceanographic settings, both display coherent orbital-scale productivity variability, indicating that basin-scale atmospheric forcing outweighed regional east–west hydrographic contrasts. In the western tropical IO, orbitally driven changes in trade winds regulated thermocline depth and nutrient supply, whereas, in the eastern tropical IO, summer monsoon intensity controlled the strength of the Java upwelling system (Andruleit et al., 2008; Tangunan et al., 2017). These results suggest that long-term productivity in both regions was primarily governed by orbitally driven atmospheric circulation despite differences in local oceanographic processes. In contrast, productivity in the NEIO was primarily regulated by freshwater-controlled stratification rather than wind-driven upwelling. Variations in precipitation, river discharge, and ITF transport governed upper-ocean stability and nutrient exchange, producing a productivity response opposite to that observed in the western and eastern tropical IO. These contrasting hydrographic controls explain the opposing PC1 loadings despite the influence of the same large-scale G–IG climate forcing.

Unlike PC1, PC2 captures secondary hydroclimatic variability superimposed on the long-term climatic signal. The PCA biplot (Fig. 7) reveals two complementary expressions of monsoon variability. PC2 is closely aligned with the Indian Monsoon Stack of Clemens and Prell (2003), which primarily reflects the atmospheric component of the ISM through wind-driven upwelling in the AS. In contrast, the Chinese speleothem composite δ¹⁸O record (Cheng et al., 2016), the δ¹⁸Osw-ivc reconstruction of Lauterbach et al. (2020), and the δ¹⁸Osw record from core BoB-24 (Liu et al., 2021) cluster together, reflecting the hydrological expression of the monsoon through variations in precipitation, freshwater discharge, and surface salinity. Together, these proxy relationships indicate that atmospheric circulation and hydrological processes exerted distinct controls on productivity across the IO.

The loading pattern of the first PCA further refines the interpretation of PC2. Among the productivity records, GeoB12613-1 exhibits the strongest loading on this component, whereas Site YDY05 contributes only weakly (Fig. 4C), indicating that PC2 primarily represents variability in the western tropical IO. Previous studies have shown that orbital-scale variations in SW monsoon winds and trade wind intensity regulated thermocline depth, upwelling, and nutrient supply in this region, thereby exerting a strong influence on BP (Tangunan et al., 2017). Accordingly, the close association between PC2 and the Clemens monsoon stack suggests that this secondary productivity mode is largely linked to wind-driven oceanographic processes, while such forcing played a comparatively smaller role in the southern BoB.

The hydrological proxies provide additional insight into productivity variability in the BoB. Their close association in the PCA reflects coherent changes in monsoon precipitation, continental runoff, and sea surface salinity, which collectively regulate upper-ocean stratification. Enhanced freshwater input strengthens the low-salinity surface layer and barrier layer, suppresses vertical mixing, and limits nutrient entrainment into the euphotic zone, whereas reduced freshwater input weakens stratification and promotes vertical nutrient transport (Sengupta et al., 2006; Thushara and Vinayachandran, 2016; Da Silva et al., 2017). The separation of these hydrological proxies from the wind-based monsoon stack therefore demonstrates that the monsoon influenced productivity through different mechanisms across the basin.

These statistical relationships are consistent with the reconstructed productivity history at Site YDY05. Although intensified SW monsoon circulation likely enhanced nutrient supply through a stronger SMC, Ekman-driven mixing, and mesoscale eddy activity (Vinayachandran et al., 2002, 2004; Prasanna Kumar et al., 2004), the concurrent increase in monsoon precipitation and freshwater discharge strengthened upper-ocean stratification, restricting nutrient entrainment into the euphotic zone. As a result, enhanced freshwater stratification largely offset the positive effects of wind-driven nutrient supply, indicating that hydrological forcing exerted the dominant control on BP variability at Site YDY05. This interpretation is consistent with the weak loading of YDY05 on PC2 and its close correspondence with the hydrological proxy records (Fig. 7).

The PCA and multiproxy results are synthesized in the conceptual model shown in Fig. 8. The model illustrates how the dominant G–IG forcing (PC1) and the secondary hydroclimatic variability (PC2) interacted to regulate BP strength in the southern BoB. Specifically, G–IG changes controlled the background nutrient state of the upper ocean, whereas monsoon-driven wind forcing and freshwater-induced stratification modulated nutrient availability and export production. Together, these processes explain the reconstructed BP variability at Site YDY05.

6 Conclusions

This study presents a high-resolution coccolithophore-based reconstruction of BP variability in the NEIO over the last ∼ 53 kyr using sediment core YDY05 from the southern BoB. The PP record exhibits pronounced G–IG variability, with minimum values (∼ 108 g C m⁻² yr⁻¹) during the LGM and maximum values (∼ 412 g C m⁻² yr⁻¹) during the EH, reflecting major changes in carbon export under evolving hydrographic conditions. Principal component analysis further demonstrates that G–IG climate variability was the dominant control on basin-scale productivity across the IO, while regional hydrographic processes governed the expression of this climatic forcing. In contrast to the wind-driven upwelling systems of the western and eastern tropical IO, productivity in the NEIO was primarily regulated by freshwater-induced upper-ocean stratification and nutrient limitation, with atmospheric and hydrological components of the ISM exerting complementary secondary influences. These findings provide new insight into the mechanisms regulating BP variability in the IO and emphasize the importance of regional hydrographic processes in shaping carbon export under changing climatic conditions.