Cyclic evolution of phytoplankton forced by changes in tropical seasonality

Although the role of Earth’s orbital variations in driving global climate cycles has long been recognized, their effect on evolution is hitherto unknown. The fossil remains of coccolithophores, a key calcifying phytoplankton group, enable a detailed assessment of the effect of cyclic orbital-scale climate changes on evolution because of their abundance in marine sediments and the preservation of their morphological adaptation to the changing environment1,2. Evolutionary genetic analyses have linked broad changes in Pleistocene fossil coccolith morphology to species radiation events3. Here, using high-resolution coccolith data, we show that during the last 2.8 million years the morphological evolution of coccolithophores was forced by Earth’s orbital eccentricity with rhythms of around 100,000 years and 405,000 years—a distinct spectral signature to that of coeval global climate cycles4. Simulations with an Earth System Model5 coupled with an ocean biogeochemical model6 show a strong eccentricity modulation of the seasonal cycle, which we suggest directly affects the diversity of ecological niches that occur over the annual cycle in the tropical ocean. Reduced seasonality in surface ocean conditions favours species with mid-size coccoliths, increasing coccolith carbonate export and burial; whereas enhanced seasonality favours a larger range of coccolith sizes and reduced carbonate export. We posit that eccentricity pacing of phytoplankton evolution contributed to the strong 405,000-year cyclicity that is seen in global carbon cycle records. .Morphometric analysis of coccolith assemblages spanning the last 2,800,000 years suggests that the evolution of coccolithophores is linked to seasonality changes, paced by Earth’s orbital eccentricity with implications for the carbon cycle.

Coccolithophores precipitate half of the biogenic CaCO 3 that is exported from the open ocean 7 and their fossil platelets (coccoliths) first appeared in sediments during the Upper Triassic, around 215 million years ago (Ma). Thereafter, coccolithophores rose to dominance 8 and became a key biological modulator of the global carbon cycle through photosynthesis and calcification 9 . In the dominant Cenozoic Noelaerhabdaceae family (including Emiliania huxleyi and Gephyrocapsa), species are defined by the morphological characteristics of their coccoliths, with size being a key criterion 10 that is related to cell size 11 . For Gephyrocapsa and Emiliania, phylogenies reconstructed from gene sequences indicate that morphology-based definitions correspond to biological species 3,12 . Within a given Noelaerhabdaceae population-which is typically dominated by one species but includes several-interspecific and intraspecific changes in coccolith length and mass occur in response to environmental parameters such as carbonate chemistry 1 and temperature 2 . Studies of coccolithophore evolution have focused on geological-timescale changes in species richness and turnover 13 , coccolith carbonate accumulation 8,14 or calcification potentially driven by carbon cycle changes 15 . In addition, climate changes induced by orbital cycles (on timescales of tens to hundreds of thousands of years) strongly influence the composition of nannofossil assemblages [16][17][18] . However, so far the effects of orbital cycles on coccolithophore evolution, coccolith morphology and carbonate production have not to our knowledge been examined simultaneously.
Here we quantify the Pleistocene history of tropical Noelaerhabdaceae evolution at high resolution (around two thousand years, kyr), using coccoliths preserved in nine well-dated sedimentary sections from the Indian and Pacific Oceans cored during International Ocean Discovery Program (IODP) and International Marine Past Global Changes Study (IMAGES) expeditions (Extended Data Table 1). We use artificial intelligence microscopy to create a biometric database of over 7 million coccoliths from more than 8,000 samples (Methods). The strong similarity of morphometric patterns observed at each site (Extended Data Fig. 1) led us to build composite frequency contour plots of coccolith size and mass, representing larger-scale evolutionary change (Fig. 1a, Methods). Patches denoting high frequency of a particular size correspond in many cases to described acmes of Noelaerhabdaceae species 19-21 or proposed evolutionary events 3 (Fig. 1). The most recent evolutive phase, which started around 550 thousand years ago (ka), is attributed to a radiation event and the emergence of new Gephyrocapsa species, on the basis of a genetic study of extant taxa and its temporal correlation to low-resolution coccolith morphometric data 3 . Over the Pleistocene, average coccolith size shows an increase that corresponds to a gradual shift in dominance from smaller Article to larger coccoliths (Fig. 1b). On orbital timescales, global ice volume and deep-sea temperature as represented by benthic foraminiferal δ 18 O show a dominance of 41-kyr and later around 100-kyr glacial-interglacial cycles 22 (Fig. 1c). By contrast, average coccolith length follows a regular cycle that is highly coherent (greater than 99.9%) with the orbital eccentricity periods of 405 kyr (e405) and of 124 and 95 kyr (e100) 23 (Extended Data Fig. 2a), with larger average size occurring at high eccentricity with a slight time lag (Fig. 1b).
Average size or mass of coccoliths in a Noelaerhabdaceae population may vary because of macro-and/or microevolution, or because ecological changes modulate the relative abundances of species in different size ranges. To build a metric that describes only species evolution, we remove the effect of relative abundance changes related to ecology 24 by formulating a morphological divergence index (MDI), calculated as the difference in average coccolith mass between two size classes-larger and smaller than 3 µm (Methods). Thus, MDI quantifies morphological divergences of species over time through evolution, and could be driven by changes in size or degree of calcification (see Fig. 2 for a conceptual explanation). Noelaerhabdaceae coccolithophores spread rapidly throughout the oceans and are often cosmopolitan, resulting in the same species being present in many regions, but with different relative abundances 12,19 . MDI varies independently of regional ecological specificities, and MDI records from sites in distinct oceanographic biomes 25 and climatic regimes (for example, warm pool, monsoon-dominated; Extended Data Table 1) are highly intercorrelated, all showing significant e405 and e100 periods (Extended Data Figs. 1, 2). Therefore, we produce a composite MDI stack, which preserves the high resolution of each dataset (Fig. 1e, Methods). The MDI stack, interpreted as reflecting evolutionary changes in morphological diversity, shows strong 405-kyr pacing throughout the Pleistocene irrespective of glacial-interglacial background state. Cross-spectral analysis indicates significant (greater than 90%) coherency between the stack and Earth's eccentricity periods since 2.8 Ma (Fig. 1d). This pattern cannot be the result of differential dissolution on coccolith morphology (Methods) and in contrast to MDI, Pleistocene deep-sea CaCO 3 dissolution generally follows glacial-interglacial cycles 26 . Similarly, coccolith morphological evolution appears not to be responding directly to physical parameters covarying with global ice volume, such as sea level or ocean temperature. Although eccentricity forcing on coccolithophore productivity has previously been suggested 27,28 , our new dataset reveals that eccentricity cycles instead forced the evolution of the Noelaerhabdaceae.
Cyclic coccolithophore evolution may have affected the ocean carbon cycle through coccolith carbonate production and burial in sediments 14,29 . Coccolithophores produce large amounts of calcite during blooms 27,30 , and sediments are often dominated by a few opportunistic species, for example E. huxleyi (0-90 ka) 19 and Gephyrocapsa caribbeanica (280-570 ka) 20 in the late Pleistocene. We estimate the mass accumulation rate of Noelaerhabdaceae coccoliths (NoMAR) in our cores and produce a stacked record (Fig. 1g, Methods). Noelaerhabdaceae coccoliths represent on average half of the total calcareous nannoplankton mass in our studied cores (Extended Data Table 2). The two components of NoMAR, coccolith flux and average mass, are separated in Extended Data Fig. 3. This reveals that NoMAR is primarily driven by changes in coccolith flux, and that flux and mass often have opposing effects on NoMAR as medium-sized, lighter species (for example, E. huxleyi and G. caribbeanica) contribute the most to coccolith carbonate export. Thus, higher NoMARs when mid-size opportunistic species dominate often correspond to lower MDI values (Fig. 1e, g). The dominance of these opportunistic species coupled with high coccolithophore accumulation in sediments during eccentricity minima is also recorded in the extra-tropics 27 . In contrast to MDI, local ecological conditions affecting productivity and export-and possibly water depth affecting coccolith accumulation-also influence NoMAR, so a linear relationship between the two is not expected. Although it is impossible to quantify the relative effects of these factors, common trends between sites emerge despite different absolute values and these are reflected in the NoMAR stack. Thus, NoMAR combines global evolutionary and local ecological drivers of calcite production, whereas MDI should exclusively record evolution. Nevertheless, the NoMAR composite record shows strong eccentricity periodicities that are significantly coherent with MDI throughout the Pleistocene (Fig. 1f), showing a strong imprint of coccolithophore morphological evolution on carbonate production and burial.

MDI and long-term seasonal variations
We hypothesize that the MDI index responds to variations in the amplitude of tropical seasonality. In low latitudes, seasonal contrast is related to the eccentricity of Earth's orbit 23,31 both directly, because the ellipticity of the orbit determines the distance between the Sun and the Earth during each season, affecting radiation intensity, and indirectly, because eccentricity modulates the effect of precession on seasonal insolation contrast. Seasonal contrast is greater during periods of high eccentricity. To our knowledge, the eccentricity-paced rhythm of surface-ocean seasonality that dominates MDI has not been documented previously because most proxies record integrated annual average conditions or a specific season. In the modern intertropical ocean, large seasonal changes in the properties of the upper water column (for example, mixed-layer depth and nutrient availability) are associated with the seasonally reversing monsoon systems and latitudinal migrations in the intertropical convergence zone. The seasonal succession of coccolithophore species, a characteristic of phytoplankton ecology, is indicative of their adaptation to the different ecological niches created by seasons 24 . In the modern ocean, the highest phytoplankton diversity is found in the tropical band, a pattern probably related to high temperatures and stable conditions, whereas seasonal species turnover is highest at mid-latitudes because of a strong seasonal temperature contrast 32 . Intra-annual dynamics of net primary production (NPP) are good descriptors of the range of oceanographic niches and biomes 25 , because NPP represents the integrated biological response to all of the changes forced by the ocean-atmosphere coupled system. To demonstrate the effect of orbital configuration on NPP seasonality and therefore niche availability, we simulated monthly oceanic NPP using the fully coupled IPSL-CM5A2 model 5 , which includes the ocean biogeochemistry model PISCES-v2 6 , for seven early Pleistocene time intervals that cover a large eccentricity spectrum with different precession conditions but with similar ice volume and obliquity (Extended Data Table 3, Fig. 1b). The results of these simulations for the tropical Indian and western Pacific Oceans show that the seasonal range of NPP increases with eccentricity, a trend that parallels the eccentricity sorted values of MDI in our Plio-Pleistocene time series (Fig. 3).
In the simulations, the increase in amplitude of the NPP seasonal cycle (Fig. 3a, Extended Data Fig. 4a-e) is primarily driven by higher productivity during boreal summer, especially in the eastern Indian ocean. This increase is forced by the modification of atmosphere-ocean dynamics in response to variations in the amplitude and seasonality of insolation forcing (Extended Data Fig. 5a-c). Eccentricity acts on sea-level pressure over continental Asia (Extended Data Fig. 5d-h) through insolation, inducing modifications of sea-level pressure gradients and low-level wind circulation over the Indo-Pacific Warm Pool (IPWP) (Extended Data Fig. 4f-h). Changes in atmospheric dynamics are responsible for regional and seasonal enhancement of NPP at high eccentricity (Extended Data Fig. 4a-c), either through the generation of anomalous upwelling along the equator (southwest of India) or the modification of the hydrological cycle that create more favourable conditions for intense vertical mixing (Extended Data Fig. 6a, c), depending on precession. Overall, those localized increases in the amplitude of the seasonal cycle lead to a less homogeneous upper ocean in the IPWP region at high eccentricity (Fig. 3a, Extended Data Fig. 4a-c). We propose that during times of high eccentricity, the higher seasonal  Article range of NPP in our model simulations (representing up to 100% of mean annual NPP) is indicative of more diverse ecological niches to which coccolithophores can adapt. A greater diversity of ecological niches when seasonality is high 25 leads to a larger number of species because Noelaerhabdaceae adaptation is characterized by the adjustment of coccolith size and degree of calcification to thrive in the new environments 1,2 .

Eccentricity lag, origination and dominance
Coccolith morphological diversity clearly responds to eccentricity (Fig. 1); however, in stark contrast to Plio-Pleistocene climate proxy records 22,33 and coccolithophore assemblage dynamics 16-18 , precession and obliquity cycles are absent from the 2-kyr resolution MDI records (Extended Data Fig. 2c-k). These cycles could have been smoothed out by the evolutionary process acting as a low-pass filter, providing an explanation for the phase lag observed between eccentricity and Noelaerhabdaceae morphology (Fig. 1b, Extended Data Fig. 2a, b). We know that speciation events spread rapidly throughout the oceans 12,19,21 , and species dominance takes longer, as exemplified by E. huxleyi. This species appeared at 290 ka but did not become dominant until 90 ka 19 , during an intense low-eccentricity interval two e100 cycles later, when it gained a competitive advantage over Gephyrocapsa oceanica and Gephyrocapsa ericsonii (Extended Data Fig. 7). The delay between species appearance and dominance could therefore be intrinsic in smoothing out variability at precession and obliquity timescales in the MDI record (Fig. 1b, Extended Data Fig. 7). At present, the lack of precise ages for the first occurrences of Noelaerhabdaceae species in the fossil record precludes us from testing this hypothesis further.
The eccentricity lags and transfer of spectral power from high to low frequencies described here are analogous to modelling results in a previous study of deep-time carbon cycle variations on orbital timescales 34 , hinting that coccolithophores may drive-rather than just respond to-carbon cycle changes.

Coccolithophores and the global carbon cycle
The persistence of e100 and e405 cycles in Cenozoic and Mesozoic records of the ocean carbon cycle (for example, per cent CaCO 3 and foraminiferal δ 13 C), independent of glacial-interglacial climate state, attests to the importance of biogeochemical processes operating at these timescales throughout Earth's history 35,36 . During the Pleistocene, Mediterranean surface δ 13 C records document e405 cycles more faithfully than do deep open-ocean records, suggesting a low-latitude climatic origin of this signal 28 . Chemical weathering has been suggested as a potential modulator of the ocean carbon cycle on 400-kyr timescales 37 . Similar to our coccolith records, a notable phase lag between δ 13 C and eccentricity is observed in the e405 band, which has been explained by the long residence time of carbon in the oceans and resultant transfer of energy from precession into eccentricity bands through a non-linear process 34,35 .Previous coccolith records spanning up to around 1 million years (Myr) have linked coccolithophore production to eccentricity forcing 17,27 . Yet changes reconstructed at our low-latitude sites cannot be explained by the hypothesis that eccentricity-driven changes in growing season length are responsible for the approximately 400-kyr cycle in coccolithophore production 27 . Our data and model results support the alternative hypothesis that changes in seasonality caused by the eccentricity of the Earth's orbit paced tropical We consider an evolutionary sequence in which species A (intermediate morphology) evolves into species B (smaller) and C (larger). In epoch 2, B and C have equivalent proportions, and in epoch 3, they fluctuate in relative abundance between 25% and 75%. The evolutionary event between epochs 1 and 2 (red shading) is not detected in mean population morphology (for example, size, mass). In epoch 3, fluctuating ecology produces population dynamics detected in mean morphology (grey shading). With a biometric boundary of 2 units, MDI jumps from 0 to 2 from epoch 1 to 2, thus it is diagnostic of an evolutionary event. In epoch 3, MDI remains stable despite fluctuating assemblage composition. In this idealized example, average population biometry is related to ecology and MDI to evolution. b, MDI calculated for IODP site U1485 (example dataset). Right, average Noelaerhabdaceae coccolith mass (smoothed using a Loess function). Middle, two size classes are created: coccoliths shorter and longer than 3 µm (grey histograms). MDI is the difference between the average log(mass) of each class (light and dark grey dots on histograms Noelaerhabdaceae evolution and production throughout the Pleistocene. Although these changes clearly affect carbonate accumulation patterns (Fig. 1g), coccolithophore productivity alone cannot be responsible for the expression of long eccentricity cycles in climate records because they are only one constituent of the phytoplankton. Other phytoplankton groups, some with little or no fossil record, may also have been similarly influenced by variations in tropical seasonality on these timescales. In this case, the effect of changes in the ratio of exported organic carbon production to carbonate mineral production, known as the rain ratio 38 , may have been strong enough to modulate the carbon cycle. The cyclic evolution of calcifying phytoplankton on eccentricity timescales in response to seasonality documented here provides evidence in support of the hypothesis that biosphere productivity must have responded to changes in solar insolation 35,37 to explain the strong e405 signature in carbon cycle records.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-021-04195-7. Miocene 386)), compared to the seasonal NPP contrast (maximum minus minimum month) from seven numerical simulations (Methods). Model points represent a regional mean of the entire map area in a. As in Fig. 1b, orange diamonds are model runs with perihelion in December (P max : 2,222 ka, 2,265 ka and 2380 ka), green stars are runs with perihelion in June (P min : 2,230 ka, 2,346 ka, 2,369 ka, and 2,395 ka)-illustrating that eccentricity has a much larger effect on seasonality than precession at a given eccentricity. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acquisition of coccolith data
Over 8,000 samples were extracted from sediment cores for coccolithophore analysis at depth intervals to achieve a high stratigraphic resolution (0.5 to 2.3 kyr; Extended Data Table 1). Samples were prepared using the settling method 39,40 : sediments were disaggregated in water and suspensions were settled onto a 12 × 12-mm cover slip and mounted with Norland Optical Adhesive 74, with 8 cover slips per microscope slide. Some samples were prepared as independent duplicates. Two slides (16 samples) were placed onto the stage of an automated polarizing microscope (Leica DM6000). After auto-focusing, 165 contiguous fields of view (with an area of 125 × 125 µm each) were imaged in each sample using a black and white SPOTFLEX camera (Diagnostic Instrument). SYRACO, a software program based on an artificial neural network 41 , identified all specimens belonging to 33 groups of coccolithophore taxa in the images 42 . The gephyrocapsid specimens-the dominant group studied here-were classified into six distinct classes that were merged into one group. On average, 888 Noelaerhabdaceae coccoliths were identified in each sample. Among other morphometric parameters, the size and mass of the coccoliths were measured. Coccolith mass is measured using birefringence, following published state-of-the-art methods 40,43 . The use of artificial intelligence in this type of work is essential because it is the only way to measure such a large number of specimens (more than 7 million) in a reasonable time, and thus obtain the high-resolution multi-site records required for this study. The pattern recognition was performed with a structured multi-layer neural network called SYRACO, written in C++ by D. Dollfus 44 . The input image of 64 × 64 pixels is connected to the output (class name) by three convolutional layers of 1,764, 360 and 80 neurons with no shared weights, which induces a long computing time. The advantages of this structure are discussed in a previous study 45 . To mimic the dynamic process of human recognition, in parallel to the second and third convolutional layers, there are 3 small neural networks of 20 neurons each, called motor layers, that perform simple image transformations from 5 possibilities: rotation, translation, symmetry, contrast and dilation. These parallel neural networks enhance the efficiency of the pattern recognition by 50% (ref. 41 ) with an accuracy above 95% (based on more than 5,000 test images). In his PhD thesis, N.B. increased the number of calcareous nannofossil species recognized by SYRACO to include most Cenozoic species and grouped them into 49 morphological classes 42 . The number of false positives (non-coccolith particles of calcareous debris such as broken foraminifera, micrite and broken coccoliths) has been reduced in SYRACO by adding a second pattern recognition level after the SYRACO artificial neural network (ANN), based on a random forest algorithm 46 . This cross-checking is more robust because it results in only 5% of false positives, compared to around 50% before 42 . In this work, we combine the 49 morphological classes into only 5 groups and work essentially with one of these, the Noelaerhabdaceae. From the confusion matrix produced by the analysis of 6,888 images (ref. 42 , table 1, p.109), the percentage of successful identifications for those 5 taxonomic groups is 96% for Noelaerhabdaceae; 91% for Coccolithales; 90% for Syracosphaearales and Zygodiscales (grouped together); and 88% for other coccolith taxa. Florisphaera profunda coccoliths are recognized at a rate of 98% (ref. 42 ). Most of the losses can be explained by the quality of the captured image, owing to some particle in a large image being out of focus or luminosity and contrast problems, or aggregation of particles. We progressively solved some of these problems by developing new optical methods 40,43 and by changing the pre-processing (for example, refining image segmentation); this increased the number of recognized coccoliths without changing the proportion of the different species. Because we were satisfied with its performance, we did not test other architectures of SYRACO such as increasing the number and the size of the convolutional layers. The goal of SYRACO was to provide a robust and rapid coccolith extractor compatible with commercial computer performance during development in the late 1990s and early 2000s. In this work, SYRACO was processed on a Dell Precision T7910 with 2 Xeon processors (2.3 GHz) of 20 cores each and 64 GB of memory, with Windows as the operating system.

Construction of composite frequency contour plot of coccolith size
Measurements were grouped into morphological bins of 0.1 µm for coccolith length in every sample. Samples were binned into 30-kyr time windows in each core, chosen such that it is larger than the length of a precession cycle (23-19 kyr). This will prevent any bias in the size and mass distribution resulting from changes in the relative abundance of large versus small gephyrocapsid species on precessional timescales 59 . Another advantage of using a 30-kyr time window is that the high number of measurements included in each bin (on average 16,650 measurements) makes it extremely precise but easier to discern trends. To standardize each time window at each site, the numbers contained in each bin are divided by the total number of coccoliths in that time window and multiplied by 100. To stack the records and produce the frequency-density plot of size (Fig. 1a), samples in each core were grouped into 30-kyr bins, standardized (%) and merged into a single stack. Frequency contour plots for size and mass (the latter not shown) show near-identical trends and variability. The distribution of coccolith mass values is skewed toward heavy values. We therefore used the logarithm of the mass to obtain a symmetrical mass distribution before binning (0.05 log(pg) bins) and stacking as for length.

Note on taxonomy of the Pleistocene Noelaerhabdaceae
The genus concept in the Pleistocene Noelaerhabdaceae is rather straightforward 10,60 : Emiliania presents 'T-shaped' elements in its distal shield, Gephyrocapsa presents a bridge in its central area, Pseudoemiliania presents slits in its distal shield and Reticulofenestra has none of these features. Two of these features may be present on Gephyrocapsa (for example, Gephyrocapsa protohuxleyi). The species concept is much more complex 10,60,61 . It essentially depends on coccolith size, ratio of central area opening to coccolith size, and bridge angle (for Gephyrocapsa). All of these features are continuous rather than discrete parameters and therefore often present a continuum between species. In Fig. 1a, b, the size density plot and the average size plot illustrate how size is a variable feature. One of the main taxonomic parameters of this group (that is, size) is constantly evolving, complicating the common use of a size-based typological species concept. A previous study 3 indicates that all extant species evolved from G. caribbeanica around 550,000 years ago, implying a rapid (less than 0.55 Myr) species turnover. The Noelaerhabdaceae family is therefore rapidly evolving genetically and morphologically. In this paper we do not intend to dispute taxonomy or species concepts. Given that it is difficult to follow genetic and typological species through time in this family, we prefer to discuss morphological evolution in a taxon-free manner. In the described taxonomies of this family 60 , there remains however a clear cut-off at around 3 µm between coccoliths of smaller and larger Noelaerhabdaceae species, which is why we choose this boundary to develop the MDI concept.

MDI
To quantitively capture the history of biological evolution within a group of species with a biometric tool, it is necessary to build a metric that it is as independent as possible from the population dynamics of the different species relative to the others. This is because the biometry of a multi-species population is greatly influenced by its population dynamics: the relative success of one species in one particular biotope will affect the average biometry of the entire population. The average biometry is therefore greatly influenced by species adaptation to biotopes and will not be necessarily diagnostic of biological evolution. A way to limit the influence of the relative abundance of species in a sample is to measure the difference between the biometric means of the considered species. A simple index can be designed to parametrize a biometric boundary between two populations: the arithmetic distance between the two population means. A schematic example of an evolutionary sequence as recorded by MDI is described in Fig. 2a.
Therefore, the MDI is not designed: (1) to trace the spatial variation of ecological parameters such as seasonality. This is because a new morphological trait will spread rapidly (on geological timescales) in the ocean if it is successful. The MDI will be distributed evenly wherever the species are present; or (2) to describe a physiological adaptation to a fluctuating environment. MDI is designed (1) to trace the morphological evolution of a small group of species; (2) to trace temporal variations of ecological parameters at the large geographical scale that can lead to the evolution of new morphological traits.

MDI designed for Plio-Pleistocene Noelaerhabdaceae
In each sample, individual coccoliths were divided into two size classes: coccoliths longer and shorter than 3 µm. The average log(mass) is calculated in both classes. MDI is the difference between the two averages (Fig. 2b). The size of 3 µm corresponds to the best cut-off value between the two modes (2.8 and 3.9 µm) of the size distribution. Other size cut-offs (2.9, 3.25 and 3.5 µm) as well as a mass cut-off (3.16 pg) were tested, without large differences in the resulting MDI values and temporal trends (Fig. 2b). The records are resampled (by linear integration) at 2-kyr intervals for further analysis (time series, statistics, and stacking). A stacked record composed of all records is calculated for each time window. This stacked MDI reflects the variability seen in all individual records (Extended Data Fig. 1). Because not all records cover the entire 2.8-Myr interval (3 records are over 2.3 Myr long, 3 are between 0.7 and 1.8 Myr long and 3 are 0.4-0.6 Myr long), the stack is composed of more records in the younger part than in the older part. Because of the phase lag between MDI and eccentricity we use band-passed eccentricity (red line in Fig. 1b) to sort and bin MDI values used for Fig. 3b. Finally, because the relative abundances of small versus larger Noelaerhabdaceae are not considered during the calculation of the MDI, any preferential dissolution of smaller more fragile coccoliths would not affect the MDI, as it represents the difference between the mean masses of the two size groups. A negligible effect of carbonate dissolution on MDI is supported by the fact that species-specific mean coccolith mass is conserved in dissolution experiments 62 , and by the similarity between MDI records regardless of core depth in the range of around 1,100-3,000 m (all well above the Pleistocene Pacific carbonate saturation horizon) (Extended Data Fig. 1).

Mass accumulation rates
Mass accumulation rates of Noelaerhabdaceae coccoliths (NoMARs) were estimated in seven cores (all cores excluding MD05-2530 and U1446), for which a quantitative sample preparation technique was applied 40 . The samples were prepared by settling onto coverslips that were weighed before and after settling, the weight difference providing the amount of sediment deposited. The number and the mass of the Noelaerhabdaceae is estimated by SYRACO. From these quantities it is possible to estimate the weight of Noelaerhabdaceae per gram of sediment. NoMAR is obtained by multiplying weight per gram by the sedimentation rate and the dry bulk density of sediment. The dry bulk density was estimated from continuous measurements of wet bulk density from gamma ray attenuation (GRA) and transformed by the linear relationship for each site between discrete shipboard measurements of wet bulk density and dry bulk density 63 . NoMARs for the 7 cores were stacked together after resampling each record at 2-kyr intervals (Fig 1g), using the same method as for MDI. Other stacking methods, such as assembling Loess-detrended records, were tried and produced consistent results. Differences exist between individual NoMAR records owing to regional difference in coccolithophore productivity, export dynamics and core depth (although only 2 cores, MD97-2140 and U1443, were retrieved from sediments deeper than 2,000 m). However, three common patterns emerge in all individual records: an increasing trend in NoMAR towards the present, a stepwise increase at around 1.1 Ma and the clear presence of eccentricity cycles. Noelaerhabdaceae coccolith flux (Extended Data Fig. 3b) is calculated as the number of coccoliths per gram of sediment multiplied by the sedimentation rate, and is the main driver of the step increase in NoMAR at around 1.1 Ma.

Time-series analysis
Time-series analyses were performed using the software packages Analyseries 64 and Acycle 65 on detrended records. Cross-spectral analyses were performed in Analyseries using Blackman-Turkey transforms 66 . For evolutive cross-spectral analyses (Fig. 1d, f) a window of 500 kyr (250 data points) and a step of 100 kyr was used. Coherence values above 0.56 are above the 80% confidence level. Spectral analyses were performed with the multi-taper method (MTM) 67 with both evolutive and entire series (Extended Data Fig. 2). Spectral properties are similar in all individual MDI records, and show that the absence of precession (19-23 kyr) and obliquity (around 41 kyr) is not a result of chronological bias in constructing the stack that would have smoothed the record. Each record has a resolution of around 2 kyr with a precise independent age model. The absence of precession and obliquity is therefore a common and robust feature of all of the MDI series as well as the stacked record.

Low-pass filters
We designed second-order low-pass filters to reproduce the effect of the time needed for a new evolved species to fully succeed (200 kyr for E. huxleyi). We transformed the following classical low-pass filter complex transfer function H: where A is the amplitude, Q is the quality factor, ω is the angular frequency 2πf (f the frequency)) in its associated differential equation: The FAD of this species is well documented because its characteristic T-shaped elements are a morphological feature that appeared suddenly, without gradation. The other Gephyrocapsa species have been described using criteria that are subject to gradation between species: coccolith length, size of the central opening, and orientation of the bridge 68 . For example, the FAD of a typical G. ericsonii, (a small gephyrocapsid) that appeared at about the same time as E. huxleyi 3 is not reported precisely because it evolved progressively from G. caribbeanica (a mid-size species). It is interesting to note that the FAD and the BA of E. huxleyi occurred similarly in times of eccentricity decrease, but two cycles apart. The intermediate cycle may have been too high to allow E. huxleyi to begin its dominance. This may not have been the case for other species under different orbital configurations. This is why we did a filter with a different configuration, which produces a delay of about one eccentricity cycle between a FAD and a BA. To express the response of those filters, we built their Bode magnitude plots, expressing the frequency response, and their Bode phase plots, expressing the phase shift (Extended Data Fig. 7b, c).

Model description
To simulate changes in NPP related to changes in eccentricity we used the Earth System Model IPSL-CM5A2 5 that simulates the interactions between ocean, atmosphere, land and ice. The following section provides a brief description of model components and experimental set-up. We then describe the model behaviour at low eccentricity and discuss how the large-scale ocean-atmosphere circulation at high eccentricity in our simulations compares to previous modelling studies.
The IPSL-CM5A2 coupled model is a combination of the LMDZ5A atmospheric model 69 72 ) and a sea-ice thermodynamics model (LIM2 73 ), as well as a biogeochemistry model (PISCES-v2 6 ), and has an horizontal resolution of 2° by 2° (refined to 0.5° in the tropics) and 31 vertical levels, the thickness of which increases from 10 m at the surface to 500 m at the bottom. The atmospheric grid has a horizontal resolution of 1.875° in latitude by 3.75° in longitude with 39 vertical levels. The ocean-atmosphere coupling is ensured by the OASIS coupler 74 that interpolates and exchanges variables between the two components. Detailed descriptions of IPSL-CM5A2 and performances in simulating pre-industrial climate can be found in previous reports 5,75 . PISCES-v2 simulates the main oceanic biogeochemical cycles (C, P, Si, N and Fe) and has a simple representation of the lower trophic levels of the marine ecosystem 6 , with two phytoplankton (nannophytoplankton and diatoms) and two zooplankton (micro-and mesozooplankton) size classes and five limiting nutrients (Fe, NO 3 − , NH 4 + , Si and PO 4 3− ). Phytoplankton growth is controlled by nutrients, light availability and water temperature. In the version of the model we used, river supply to the ocean of all elements apart from dissolved inorganic carbon and alkalinity is taken from the GLOBAL-NEWS2 datasets 76 and does not vary from one simulation to another. Model parameterizations are detailed in a previous report 6 .
Simulations were performed for seven early Pleistocene time slices and differ only by their respective orbital parameters (Extended Data Table 3, Fig. 1b). The time slices were chosen to target the signal produced by the 405-kyr eccentricity cycle. Land-sea mask, ice-sheet configuration as well as CO 2 and other greenhouse gases concentrations are set at pre-industrial values. Each simulation was started from the same equilibrated pre-industrial simulation 5 and was run for 500 years. NPP is integrated over the whole water column and averaged over the last 100 model years.
At low eccentricity (E min P max and E min P min ) the eastern Indian ocean surface dynamics are forced by the summer westerlies that blow northward over the Bay of Bengal (Extended Data Fig. 4f), associated with high precipitation over India and the Himalayan foreland region, whereas strong easterlies are recorded south of the equator. Winds force strong westward surface currents along the equator and south of Sumatra that generate upwelling (Extended Data Fig. 7b). The latter advects nutrients to the surface (Extended Data Fig. 7a) and triggers high productivity during summer. This peak productivity contributes to the strong seasonal cycle in this region. The winds reverse during boreal winter, triggering a second productivity bloom of lesser intensity (not shown). The productivity minimum is recorded during late spring when low-level winds along the equator are weak westerlies that favour downwelling and prevent strong convective mixing, which results in lower nutrient content within the surface layer of the ocean. The seasonal cycle of productivity in this region is very similar to the cycle simulated for the present-day equatorial Indian Ocean 5 .
During high eccentricity periods at precession minima (maxima), increasing (decreasing) boreal summer insolation (Extended Data  Fig. 5b, c) is responsible for increasing (decreasing) sea-level pressure over continental Asia (Extended Data Fig. 5d-h). Induced modifications of sea-level pressure gradients over the tropical Indian and Pacific Oceans in turn translate into changes in the low-level wind circulation over the IPWP (Extended Data Fig. 4f-h). Anomalous easterlies at precession minima (westerlies at precession maxima) in the equatorial region generate anomalous upwelling along the equator (southwest of India) that are responsible for the increasing nutrient content at the surface triggering large enhancement of productivity (Extended Data Fig. 6a-c). NPP is, in addition, amplified by modifications of the hydrological cycle that create more favourable conditions related to changes in salinity, water temperature and/or amount of solar radiation at the surface. At maximum precession and eccentricity, for example, higher sub-surface salinity (+0.5 to 1.6 psu) and lower temperatures (−1.2 to −2 °C) in the western Bay of Bengal (Extended Data Fig. 7c) reduce stratification of the upper-water column, which favours vertical mixing and contributes to enhanced productivity. The simulated patterns of atmosphere-ocean circulation and surface ocean physical state (Extended Data Figs. 4f-h, 6b, c) are in line with previous modelling studies under similar orbital configurations 77-79 . In addition, our simulations also illustrate how these changes affect the seasonal productivity cycle. The increasing amplitude of the seasonal cycle in the surface ocean at high eccentricity is probably not limited to the IPWP area. For example, another study 80 also simulates enhancement of the surface ocean temperature cycle at high eccentricity in the Eastern Equatorial Pacific, with higher amplitude than in the Western Equatorial Pacific.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

Data availability
All coccolith morphological data, as well as all model outputs described in the paper (including NPP and main oceanic and atmospheric variables) are archived at the SEANOE open access data repository: https:// doi.org/10.17882/84031. LMDZ, XIOS, NEMO and ORCHIDEE are released under the terms of the CeCILL license. OASIS-MCT is released under the terms of the Lesser GNU General Public License (LGPL). IPSL-CM5A2 source code is publicly available through svn, with the following commands line : svn co http://forge.ipsl.jussieu.fr/igcmg/ svn/modipsl/branches/publications/IPSLCM5A2.1_11192019 modipsl ; cd modipsl/util ; ./model IPSLCM5A2.1. The mod.def file provides information regarding the different revisions used, namely: NEMOGCM branch nemo_v3_6_STABLE revision 6665; XIOS2 branchs/xios-2.5 revision 1763; IOIPSL/src svn tags/v2_2_2; LMDZ5 branches/IPSLCM5A2.1 rev 3591; branches/publications/ORCHIDEE_IPSLCM5A2.1.r5307 rev 6336; and OASIS3-MCT 2.0_branch (rev 4775 IPSL server). The login/ password combination requested at first use to download the ORCHI-DEE component is anonymous/anonymous. We recommend that you refer to the project website: http://forge.ipsl.jussieu.fr/igcmg_doc/ wiki/Doc/Config/IPSLCM5A2 for a proper installation and compilation of the environment.