Luc Beaufort
email: beaufort@cerege.fr

Clara T Bolton

Anta-Clarisse Sarr

Baptiste Suchéras-Marx

Yair Rosenthal

Yannick Donnadieu

Nicolas Barbarin

Samantha Bova

Pauline Cornuault

Yves Gally

Emmeline Gray

Jean-Charles Mazur

Martin Tetard

Cyclic evolution of phytoplankton forced by changes in tropical seasonality

The effect of global climate cycles driven by Earth's orbital variations on evolution is poorly understood because of difficulties achieving sufficiently-resolved records of past evolution. The fossil remains of coccolithophores, a key calcifying phytoplankton group, enable an exceptional assessment of the impact of cyclic orbital-scale climate change on evolution because of their abundance in marine sediments, and because coccolithophores demonstrate extreme morphological plasticity in response to the changing environment 1,2 . Recently, evolutionary genetic analyses linked broad changes in Pleistocene fossil coccolith morphology to species radiation events 3 . Here, using high-resolution coccolith data, we show that during the last 2.8 million years coccolithophore evolution was forced by Earth's orbital eccentricity with rhythms of ~100,000 years and 405,000 years -a distinct spectral signature to that of coeval global climate cycles 4 .

Simulations with an Earth System Model 5 including the marine carbon cycle 6 demonstrate that eccentricity directly impacts the diversity of ecological niches occurring over the annual cycle in the tropical ocean. Reduced seasonality favours species with mid-size coccoliths that bloom year-round, increasing coccolith carbonate export and burial. We posit that eccentricity pacing of phytoplankton evolution contributed to the strong 405,000-year pacing seen in records of the global carbon cycle.

 , or calcification potentially driven by carbon cycle changes 13 . However, to date a lack of records that are both long-term and high-resolution has precluded an understanding of the impact of orbital cycles (on timescales of tens to hundreds of thousands of years) on coccolithophore evolution and carbonate production.

Here, we quantify the history of tropical Noëlaerhabdaceae evolution throughout the Pleistocene at high resolution (~2 thousand years, kyr), using coccoliths preserved in nine well-dated sedimentary sections from the tropical Indian and Pacific Oceans (Fig. 1; Extended Data Table 1; Methods). We use artificial intelligence microscopy to create an unprecedented biometric database of over 7 million coccoliths from >8000 samples (see Methods). The remarkable similarity of morphometric patterns observed at each site has led us to build evolutionary histograms of coccolith size and mass. To create these, we binned data into 30-kyr intervals to erase any precession-band variability driven by ecology 14 and to increase the statistical robustness of each point, now composed of measurements of tens of thousands of individual coccoliths (Fig. 2a; Methods). Evolutionary histograms for size and mass show identical trends and variability, we therefore only discuss the size record. The data show multiple phases of contraction alternating with phases of dilation characterised by bimodal or broader morphological distributions (Fig. 2a). Patches of high abundance of a particular morphotype correspond in many cases to described acmes of Noëlaerhabdaceae species [15][16][17] . The most recent dilation phase that started around 450 thousand years ago (ka) has been attributed to a radiation event and the emergence of new Gephyrocapsa species, based on a consensus phylogeny reconstructed from genome sequence data and its temporal correlation to low-resolution coccolith morphometric data [START_REF] Bendif | Repeated species radiations in the recent evolution of the key marine phytoplankton lineage Gephyrocapsa[END_REF] . Over the entire record, average coccolith size shows an increasing trend that corresponds to a gradual shift in dominance from smaller to larger coccoliths during the Pleistocene (Fig. 2b). On orbital time scales, coccolith size and global ice volume and temperature cycles as represented by the δ 18 O of benthic foraminifers 18 show very different temporal evolution (Fig. 2c). Instead the size contraction and dilation alternations follow a regular cycle that is highly coherent (>99.9%) with the orbital eccentricity periods of 405 kyr (e405) and of 124 and 95 (e100) 19 (Extended Data Figure 1a), with contractions (dilations) occurring at low (high) eccentricity with a slight time lag (Fig. 2b).

The average size or mass of coccoliths in a Noëlaerhabdaceae population may vary over time because of macro-and/or micro-evolution, or because ecological changes modulate the relative abundances of species in different size ranges. To remove the effect of relative abundance changes, we create a novel Morphological Divergence Index, MDI, that represents the difference in average coccolith mass between two size classes, larger and smaller than 3µm (Extended Data Figure 2; Methods). MDI records from individual sites are highly intercorrelated (Extended Data Figure 3) and all show significant e405 and e100 periods (Extended Data Figure 1), despite coming from distinct oceanographic biomes 20 and climatic regimes (e.g. Warm Pool, monsoon-dominated; see Extended Data Table 1). We therefore produce a composite MDI stack, preserving the high resolution of each dataset (Fig. 3c; Methods). Evolutionary cross-spectral analysis indicates that the stack has highly significant (>90%) coherency with the frequencies of Earth's eccentricity since 2.8 Ma (Fig. 3b,b'). This pattern cannot be the result of differential dissolution on coccolith morphology, because MDI values are similar regardless of core depths, and because in contrast to MDI, Pleistocene deep-sea CaCO3 dissolution generally follows a glacial-interglacial (G-IG) pattern [START_REF] Sexton | Onset of 'Pacific-style' deep-sea sedimentary carbonate cycles at the mid-Pleistocene transition[END_REF] . The MDI stack, interpreted as primarily reflecting evolutionary changes in morphology, shows strong 405-kyr pacing throughout the Pleistocene, irrespective of G-IG background state.

While eccentricity forcing on coccolithophore productivity has previously been suggested [START_REF] Rickaby | Coccolith chemistry reveals secular variations in the global ocean carbon cycle? Earth Planet[END_REF][START_REF] Wang | Obscuring of long eccentricity cyclicity in Pleistocene oceanic carbon isotope records[END_REF] , our new high-resolution multi-million-year dataset reveals that eccentricity cycles imprinted the evolution of the Noëlaerhabdaceae. Such cyclicity may have had important consequences for the oceanic carbon cycle on multiple time scales via coccolith carbonate production and burial in sediments [START_REF] Si | Reduced continental weathering and marine calcification linked to late Neogene decline in atmospheric CO2[END_REF][START_REF] Suchéras-Marx | Downsizing the pelagic carbonate factory: Impacts of calcareous nannoplankton evolution on carbonate burial over the past 17 million years[END_REF] . Coccolithophores are known to produce high amounts of calcite during blooms [START_REF] Rickaby | Coccolith chemistry reveals secular variations in the global ocean carbon cycle? Earth Planet[END_REF][START_REF] Holligan | A biogeochemical study of the coccolithophore, Emiliania huxleyi, in the North Atlantic[END_REF] , and sediments are often dominated by few species, for example Emiliania huxleyi (0-90 ka) 15 and Gephyrocapsa caribbeanica (280-570 ka) 16 in the late Pleistocene. We estimate the CaCO3 mass accumulation rate (MAR) of Noëlaerhabdaceae coccoliths in our cores and produce a stacked CaCO3 MAR record (Fig. 3e, Methods). This reveals that higher MARs correspond to lower MDI values when mid-size bloom species, which contribute the most to coccolith carbonate export, dominate. The dominance of these bloom species coupled with high coccolithophore production and accumulation in sediments during eccentricity minima is also recorded in the extra-tropics [START_REF] Rickaby | Coccolith chemistry reveals secular variations in the global ocean carbon cycle? Earth Planet[END_REF] , attesting to the global significance of these observations. In our records, the two main components, production of coccoliths and individual mass, often have opposing effects on MAR since medium-size and low-mass species, such as E. huxleyi and G. carribeanica, are the most productive (Extended Data Figure 4). In addition, local ecological conditions affecting productivity and export, and potentially water depth affecting coccolith preservation, likely also influence MAR, although it is not possible to quantify the relative effects of these factors. Nevertheless, the composite MAR record shows strong eccentricity periodicities that are significantly coherent with MDI throughout the Pleistocene (Fig. 3d,d'), demonstrating a strong imprint of coccolithophore evolution on carbonate production and burial.

MDI as a recorder of seasonality

We hypothesise that our new coccolith index reflects variations in the amplitude of tropical seasonality. In low latitudes, seasonal contrast is related to the eccentricity of Earth's orbit, (1) directly because the ellipticity of the orbit determines the distance between the Sun and the Earth during each season, affecting radiation intensity, and (2) 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. The seasonal succession of coccolithophore species, a characteristic of phytoplankton ecology, is indicative of their adaptation to the different ecological niches created by seasons [START_REF] Okada | Seasonal distribution of modern Coccolithophores in the western North Atlantic Ocean[END_REF] . Those seasonal niches are more diverse in a climate with stronger seasonality 20 . In the modern ocean, the highest phytoplankton diversity is found in the tropical band, a pattern probably related to high temperatures and temporal stability, whereas the seasonal species turnover is highest at mid-latitudes because of a strong seasonal contrast [START_REF] Righetti | Global pattern of phytoplankton diversity driven by temperature and environmental variability[END_REF] . Annual records of Net Primary Production (NPP) are good descriptors of the range of oceanographic niches and biomes 20 . To demonstrate the effect of orbital configuration on NPP and therefore niche availability, we simulated monthly oceanic NPP using the fully-coupled IPSL-CM5A-2 model 6 , which includes the ocean biogeochemistry model PISCES-v2 5 , for seven early Pleistocene time intervals chosen to cover a large eccentricity spectrum with different precession conditions but with similar ice volume and obliquity (lines in Fig. 2; Extended Data Table 2). The results of these simulations for the tropical eastern 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 (Fig. 4).

In the simulations, the increase in the amplitude of the NPP seasonal cycle (Fig 1; Extended Data Figure 5b-e) is primarily driven by higher productivity during boreal summer, especially in the eastern Indian ocean. This increase is forced by modification of atmosphere-ocean dynamics in response to variations in the amplitude and seasonality of insolation forcing that is amplified at extreme eccentricity values (Extended Data Figure 5g-j and6a-c). During high eccentricity periods at precession minima (maxima), increasing (decreasing) boreal summer insolation (Extended Data Figure 6b,c) is responsible for increasing (decreasing) sea-level pressure over continental Asia (Extended Data Figure 6e-h). Induced modifications of sea-level pressure gradients over the tropical Indian and Pacific Oceans in turn translate into changes in the lowlevel wind circulation over the Indo-Pacific Warm Pool, IPWP, (Extended Data Figure 5g-j). Anomalous easterlies at precession minima (westerlies at precession maxima) in the equatorial region generate anomalous upwelling along the equator (SW of India) that are responsible for the increasing nutrient content at the surface triggering large enhancement of productivity (Extended Data Figure 7a-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 relative to in other simulations (Extended Data Figure 7c) reduce stratification of the upper-water column, which favours vertical mixing and contributes to enhanced productivity. 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. 1 and Extended Data Figure 5a-c). We propose that during high eccentricity times, the higher seasonal range of NPP in our model simulations is indicative of more diverse ecological niches that coccolithophores can adapt to. A greater diversity of ecological niches when seasonality is high 20 leads to a larger range of morphotypes because Noëlaerhabdaceae adaptation is characterised by the adjustment of coccolith size and degree of calcification (therefore mass) to thrive in the new environments 1,2 .

Eccentricity lags, species origination, and dominance

The MDI reflects periods of dominance of a particular species or morphotype, rather than species origination events, and is clearly modulated by eccentricity. In stark contrast to most Pleistocene climate proxy records (e.g. 14,18 ), precession and obliquity cycles are absent from the highresolution MDI records (Extended Data Figure 3), suggesting that evolutionary processes smooth the records. Significant lags exist between eccentricity forcing and coccolith morphology (Fig. 2b), which we hypothesize result from the time taken for favourable conditions to reoccur that allow a species to rise to dominance following its first emergence. For example, E. huxleyi first appeared at 290 ka, but did not become dominant until 90 ka 15 . Yet speciation events spread rapidly throughout the oceans 15,17 , a phenomenon reflected in the small genetic diversity amongst the global Noëlaerhabdaceae population [START_REF] Filatov | Extreme Lewontin's paradox in ubiquitous marine phytoplankton species[END_REF] . This time lag could result from interspecific competition. When E. huxleyi originated, Gephyrocapsa oceanica and Gephyrocapsa ericsonii already dominated coccolithophore assemblages, and the time required for E. huxleyi to become widespread while favourable conditions persisted was not sufficient for it to become dominant. It was not until two e100 cycles later, under similar eccentricity conditions, that E. huxleyi gained a competitive advantage and rose to dominance. The delay between species appearance and dominance is thus intrinsic in the MDI record, acting as a low pass-filter with a cut-off frequency of 1/200 kyr-1 (Fig. 2b, Extended Data Figure 8). For other species, the lag might have been one or two e100 cycle, providing an explanation for the non-stationary nature of the phase lag between e100 and MDI.

Long-eccentricity forcing, coccolithophores, and the global carbon cycle

The persistence of the long eccentricity cycle e405 in Cenozoic and Mesozoic records of the ocean carbon cycle (e.g., weight percent CaCO3 and foraminiferal δ 13 C) [START_REF] Pälike | The Heartbeat of the Oligocene Climate System[END_REF][START_REF] Herbert | A long marine history of carbon cycle modulation by orbital-climatic changes[END_REF][START_REF] Kocken | The 405 kyr and 2.4 Myr eccentricity components in Cenozoic carbon isotope records[END_REF] , independent of G-IG climate state, testifies to the importance of biogeochemical processes operating at this timescale throughout Earth's history. In the Pleistocene, Mediterranean surface δ 13 C records record e405 cycles more faithfully than deep open-ocean sites, suggesting a low-latitude climatic origin of this signal [START_REF] Wang | Obscuring of long eccentricity cyclicity in Pleistocene oceanic carbon isotope records[END_REF] . Similar to our coccolith records, a significant 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 via a non-linear process [START_REF] Pälike | The Heartbeat of the Oligocene Climate System[END_REF] . Previous studies spanning up to ~1 Myr have linked coccolithophore production to eccentricity forcing 14,[START_REF] Rickaby | Coccolith chemistry reveals secular variations in the global ocean carbon cycle? Earth Planet[END_REF] . However, the hypothesis proposed by Rickaby et al. [START_REF] Rickaby | Coccolith chemistry reveals secular variations in the global ocean carbon cycle? Earth Planet[END_REF] that eccentricity-driven changes in growing season length may be responsible for the ~400-kyr cycle in coccolithophore production cannot adequately explain changes reconstructed at our low latitude sites.

Instead, our data and model results lead us to the alternative hypothesis that changes in seasonality caused by the eccentricity of the Earth's orbit paced tropical Noëlaerhabdaceae production and evolution throughout the Pleistocene. Although these changes clearly impact carbonate accumulation patterns (Fig. 3), coccolithophore productivity alone cannot be responsible for the expression of long eccentricity cycles in climate records. Other phytoplankton groups, some with little or no fossil record, must also have been influenced by variations in tropical seasonality on these timescales. In this case, the impact of changes in the ratio of exported organic carbon production to carbonate mineral production, known as the rain ratio [START_REF] Archer | Effect of deep-sea sedimentary calcite preservation on atmospheric CO2 concentration[END_REF] , may have been strong enough to modulate the carbon cycle. Our data, showing cyclic evolution of calcifying phytoplankton on eccentricity timescales in response to seasonality, provide a new perspective in support of the hypothesis of Pälike et al. [START_REF] Pälike | The Heartbeat of the Oligocene Climate System[END_REF] that biosphere productivity must have responded directly to changes in solar insolation in order to explain the strong e405 signature in carbon cycle records.

Main references

Methods

Coccolith data acquisition

Over eight thousand 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 32,33 : sediments were disaggregated in water and suspensions were settled onto a 12x12 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 at the same time onto the stage of an automated polarizing microscope (Leica DM6000). Following auto-focusing, 165 contiguous fields of view (with an area of 125 x 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 [START_REF] Beaufort | Automatic recognition of coccolith by dynamical neural network[END_REF] , identified all specimens belonging to 33 groups of coccolithophore taxa in the images [START_REF] Barbarin | La reconnaissance automatisée des nannofossiles calcaires du Cénozoique[END_REF] . The gephyrocapsid specimens, the dominant group studied here, were classified into six distinct classes that were merged into one group. On average, 888

Noëlaerhabdaceae coccoliths were identified in each sample. Among other morphometric parameters, size and mass of the coccoliths were measured. Coccolith mass is measured using birefringence, following published state-of-the-art methods [START_REF] Beaufort | Optical measurements to determine the thickness of calcite crystals and the mass of thin carbonate particles such as coccoliths[END_REF][START_REF] Beaufort | Weight estimates of coccoliths using the optical properties (birefringence) of calcite[END_REF] . 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 (>7 million) in a reasonable time, and thus obtain the high-resolution multi-site records required for this study.

Site-specific chronologies

IMAGES core MD97-2140 [START_REF] De Garidel-Thoron | Stable sea surface temperatures in the western Pacific warm pool over the past 1.75 million years[END_REF] : The age model for Site MD97-2140 [START_REF] De Garidel-Thoron | Stable sea surface temperatures in the western Pacific warm pool over the past 1.75 million years[END_REF] , on the Eauripik Rise in the Western Pacific Warm Pool, is based on tuning of a high-resolution planktonic foraminiferal Globigeronoides ruber δ 18 O record to the astronomically calibrated ODP Site 677 δ 18 O G. ruber record [START_REF] Shackleton | An alternative astronomical calibration of the lower Pleistocene timescale based on ODP Site 677[END_REF] , located in the eastern equatorial Pacific. This age model yields a chronology consistent with major micropaleontological (disappearance of G. ruber pink variety) and palaeomagnetic (Brunhes-Matuyama boundary) events [START_REF] Carcaillet | Geomagnetic moment instability between 0.6 and 1.3 Ma from cosmonuclide evidence[END_REF] .

IMAGES core MD05-2920 [START_REF] Beaufort | Shipboard Scientific Party[END_REF] :The age model for Site MD05-2920 [START_REF] Tachikawa | The precession phase of hydrological variability in the Western Pacific Warm Pool during the past 400 ka[END_REF] , on the southern bank of Manus Island, is based on ten Accelerator Mass Spectrometer (AMS) 14 C dates obtained from the surface-dwelling planktonic foraminifera G. ruber (white), and a correlation between the benthic foraminiferal δ 18 O record and the reference LR04 benthic foraminiferal δ 18 O stack 18 .

IMAGES core MD05-2930 [START_REF] Beaufort | Shipboard Scientific Party[END_REF]

Construction of composite histograms of coccolith size and mass

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. Measurements were grouped into morphological bins of 0.1 µm for coccolith length and 0.05 log10(pg) for coccolith mass 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/mass distribution resulting from changes in the relative abundance of large versus small gephyrocapsid species on precessional timescales [START_REF] Beaufort | Reconstructing Ocean History: a window into the future[END_REF] . Another advantage of using a 30 kyr time window is that the high number of measurements included in each 30 kyr bin (on average 16650 measurements) make it extremely precise but easier to discern trends. To standardise 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 evolutionary size histograms (Fig. 2a) and mass, all records are summed by bin and divided by the number of records.

Morphological Divergence Index (MDI)

In each sample, individual coccoliths were divided into two size classes: coccoliths longer and shorter than 3 µm. The average log10(mass) is calculated in both classes. MDI is the difference between the two averages (Extended Data Figure 2). 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 were tested (2.9, 3.25 and 3.5 µm) without notable differences in the resulting MDI values and temporal trends. 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 Figure 4). 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. 2b) to sort and bin DMI values used for Fig. 4.

Mass accumulation rates

Mass accumulation rates of Noëlaerhabdaceae coccoliths were estimated in seven cores (all cores excluding MD05-2530 and U1446) for which a quantitative sample preparation technique was applied [START_REF] Beaufort | Optical measurements to determine the thickness of calcite crystals and the mass of thin carbonate particles such as coccoliths[END_REF] : 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 Noëlaerhabdaceae is estimated by SYRACO. From these quantities it is possible to estimate the weight of Noëlaerhabdaceae per gram of sediment. Noëlaerhabdaceae CaCO3 MAR (only called MAR hereafter) 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 (GRAPE) and transformed by the linear relationship for each site between discrete shipboard measurements of wet bulk density and dry bulk density [START_REF][END_REF] . MARs for the 7 cores

were stacked together after resampling each record at 2 kyr intervals (Fig. 3e), 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 MAR records due to regional difference in coccolithophore productivity, export dynamics, and preservation (although only 2 cores, MD97-2140 and U1443, were retrieved from sediments deeper than 2000 m). However, three common patterns emerge in all individual records: an increasing trend in MAR towards the present, a sharp stepwise increase at ~1 Ma, and the clear presence of eccentricity cycles. Noëlaerhabdaceae coccolith flux (Extended Data Figure 5c) is calculated as the number of coccoliths per gram of sediment multiplied by the sedimentation rate.

Time series analysis

Time series analyses were performed using the software packages Analyseries [START_REF] Paillard | Macintosh program performs time-series analysis[END_REF] and Acycle 51 on detrended records. Cross-spectral analyses were performed in Analyseries using Blackman-Turkey transforms [START_REF] Blackman | The measurement of power spectra from the point of view of communications engineering-Part I[END_REF] . For evolutive analyses (Fig. 3) 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 (Fig. 2b). Spectral analyses were performed with the Multi Taper Method (MTM) [START_REF] Thomson | Spectrum estimation and harmonic analysis[END_REF] with both evolutive and entire series (Extended Data Figure 1). Spectral properties are similar in all individual MDI records, and show that the absence of precession (19-23 kyr) and obliquity (~41 kyr) is not a result of chronological bias in constructing the stack that would have smoothed the record. Each record has a ~2 kyr resolution 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 2 nd order low-pass filters in order 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 the frequency)) in its associated differential equation: The FAD of this species is well documented because its characteristic T shape 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 [START_REF] Bollmann | Morphology and biogeography of the Gephyrocapsa coccoliths in Holocene sediments[END_REF] . For example, the FAD of a typical G. ercisonii, (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. In order 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 Figure 8b,c).

Model Description

To simulate changes in Net Primary Production (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 experiments setup. We then describe the model behaviour at low eccentricity and discuss how the large-scale ocean-atmosphere circulation at high eccentricity in our simulations compare to previous modelling studies.

The IPSL-CM5A2 coupled model is a combination of the LMDz5A atmospheric model [START_REF] Hourdin | Impact of the LMDZ atmospheric grid configuration on the climate and sensitivity of the IPSL-CM5A coupled model[END_REF] , the ORCHIDEE 56 land surface model and the NEMOv3.6 oceanic model [START_REF] Madec | NEMO ocean engine. Note du Pôle de modélisation de l[END_REF] . The NEMO model includes an ocean dynamic component (OPA [START_REF] Madec | NEMO reference manual, ocean dynamics component: NEMO-OPA. Preliminary version. Note du Pole de modélisation[END_REF] ), a sea-ice thermodynamics model (LIM2 [START_REF] Fichefet | Sensitivity of a global sea ice model to the treatment of ice thermodynamics and dynamics[END_REF] ) as well as a biogeochemistry model (PISCES-v2 [START_REF] Aumont | PISCES-v2: An ocean biogeochemical model for carbon and ecosystem studies[END_REF] ) and has an horizontal resolution of 2° by 2° (refined to 0.5° in the tropics) and 31 vertical levels, whose thickness 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 [START_REF] Valcke | Coupling technologies for earth system modelling[END_REF] that interpolates and exchange variables between the two components. Detailed description of IPSL-CM5A2 and performances in simulating pre-industrial climate can be found in references 61,62. 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 [START_REF] Aumont | PISCES-v2: An ocean biogeochemical model for carbon and ecosystem studies[END_REF] , with two phytoplankton (nannophytoplankton and diatoms) and two zooplankton (micro-and meso-zooplankton) size classes and five limiting nutrients (Fe, NO3 -, NH4 + , Si and PO4 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 DIC and alkalinity is taken from the GLOBAL-NEWS2 data sets [START_REF] Mayorga | Global nutrient export from WaterSheds 2 (NEWS 2): model development and implementation[END_REF] and does not vary from one simulation to another. Model parameterizations are detailed in reference 6.

Simulations were performed for seven early Pleistocene time slices and differ only by their respective orbital parameters (Extended Data Table 2, Figure 2). The time slices were chosen in order to target the signal produced by the 405-kyr eccentricity cycle. Land-sea mask, ice-sheets configuration as well as CO2 and other greenhouse gases concentrations are set at pre-industrial values. Each simulation was started from the same equilibrated pre-industrial simulation [START_REF] Sepulchre | An Earth System Model designed for long simulation of past and future climates[END_REF] 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 (EminPmax and EminPmin) the eastern Indian ocean surface dynamics is forced by the summer westerlies that blow northward over the Bay of Bengal (Extended Data Figure 6f), associated with high precipitation over India and the Himalayan foreland region, while strong easterlies are recorded south of the equator. Winds force strong westward surface currents along the equator and south of Sumatra Island that generate upwelling (Extended Data Figure 7b). The latter advects nutrients to the surface (Extended Data Figure 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 EIO [START_REF] Koné | Physical and biogeochemical controls of the phytoplankton seasonal cycle in the Indian Ocean: A modeling study[END_REF] .

At high eccentricity, the simulated patterns of atmosphere-ocean circulation and surface ocean physical state (Extended Data Figure 6g-j and 8b, c) are in line with previous modelling study under similar orbital configurations [START_REF] Bosmans | Response of the Asian summer monsoons to idealized precession and obliquity forcing in a set of GCMs[END_REF][START_REF] Braconnot | Impact of precession on monsoon characteristics from coupled ocean atmosphere experiments: changes in Indian monsoon and Indian ocean climatology[END_REF][START_REF] Prescott | Indian monsoon variability in response to orbital forcing during the late Pliocene[END_REF] . In addition, our simulations also illustrate how these changes impact 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, Erb et al. 68 and (e) represent the anomaly of yearly maximum contrast in NPP. At EmaxPmax, the eastern equatorial Indian Ocean exhibits moderate seasonality (a) due to inhibition of the summer productivity induced by lower nutrient concentrations in this area (Extended Data Figure 8a). In this case, high productivity areas during boreal summer are shifted to south-west of India. Bottom: Late summer (JASO) low-level winds for (f) EminPmin (g) EmaxPmin (h) EmaxPmax simulations. (i) and (j) represent the anomaly in late summer low-level winds. At EmaxPmax the north-equatorial westerlies (c, e) are confined to south of 10°N due to the extension above India of the low-pressure area (Figure 1).

  Bollmann, J., Baumann, K. H. & Thierstein, H. R. Global dominance of Gephyrocapsa coccoliths in the Late Pleistocene: Selective dissolution, evolution, or global environmental change? Paleoceanography 13, 517-529, doi:10.1029/98PA00610 (1998).

FiguresFigure 1 :Figure 2 :

 12 Figures

Figure 3 :Figure 4 :

 34 Figure 3: Eccentricity, coccolith morphological divergence, and Noëlaerhabdaceae CaCO3 mass accumulation rates. a: Earth's

9 :

 9 𝑥 where y is the output series and x the input series.We solved it numerically as follow: 𝜔 9 )𝑦 5.( -𝑦 5.: )𝐴𝜔 9 : 𝑥 5Two configurations have been chosen: one produces a delay of 200 kyr for the new species' success, the second produces a delay of 130 kyr (Extended Data Figure8a). To obtain those delays, we used the following values:First case (delay of 200 kyr), 𝜔 9 = 0.021, Q = 1, A = 1Second case (delay of 130 kyr), 𝜔 9 = 0.035, Q = 1.02, A = 0.65The delay of 200 kyr corresponds to the time between the First Appearance Datum (FAD) and the beginning of the acme (BA) of E. huxleyi 17 .

Extended Data Figure 1

 1 : a: Cross-spectral analysis between eccentricity and stacked average coccolith length. Top: coherency; Bottom: phase (radian). b to j: MTM and evolutive spectral analyses (see Methods) of detrended individual MDI series resampled at 2 kyr intervals (shown on left of each evolutive analysis. Primary orbital periods are shown by red lines.

Extended Data Figure 3 :Extended Data Figure 4 :Extended Data Figure 5 :

 345 Individual MDI records for each site (black lines and points) plotted with the stacked MDI record (red line). Pearson correlation coefficients between individual sites and the stack vary between 0.71 and 0.93 (p-values are all <0.00001). Figure comparing stacked records, binned into 10-kyr intervals, of Noëlaerhabdaceae average coccolith mass (a), Noëlaerhabdaceae CaCO3 mass accumulation rates (b), and Noëlaerhabdaceae coccolith flux to the seafloor (c) (see Methods). Grey shaded areas represent four described acmes of mid-size Noëlaerhabdaceae species 15-17,70,71 . Top: Yearly maximum contrast in NPP (gC.m -2 .day -1 ) for (a) EminPmin, (b) EmaxPmin and (c) EmaxPmax. Low eccentricity values minimise the amplitude of precession variability (Extended Data Figure 7a-c), thus we only show results for minimum precession value at minimum eccentricity (EminPmin) but the reader can consider those results to be similar for the EminPmax simulation. (d)

Extended Data Figure 6 :∆)

 6 Seasonal latitudinal variations of solar radiation at the top of the atmosphere derived from the model (W.m -2 ). (a) EminPmin, (b) EmaxPmin, (c) EmaxPmax. See Extended Data Table 2 for details of orbital configurations of each simulation. Late summer (JASO) low-level winds for (d) EminPmin (e) EmaxPmin (f) EmaxPmax simulations and anomaly in late summer low-level winds (g) EmaxPmin minus EminPmin (h) EmaxPmax minus EminPmax. Low-level winds (m.s -1 ) Low-level winds (m.s -1 Extended Data Figure 7: (a) NO3 concentrations in the surface layer (0-100m), (b) upwelling velocity (averaged between 40 and 80m), (c) SST. All variables are averaged over JASO. Left : Emin, Middle : EmaxPmin minus EminPmin, Right : EmaxPmax minus EminPmax. Extended Data Figure 8: a and b: Bode plots of the 130 kyr delay filter (red) and the 200 kyr delay filter (blue) for magnitude (a) and phase (b)(see Methods). Earth's primary orbital periods are indicated by shading. c: Low-pass filter design, for the delay between First Appearance Datum (FAD) and the beginning of the acme (BA) for E. huxleyi (blue line) and other gephyrocapsids (red line). The stepped green line represents E. huxleyi's existence (filtered such that 0 means absence and 1 presence). The blue and red curves are the output series of the 2 lowpass filters described in the methods. The black curve represents coeval eccentricity values.

  ). Blue dots are model runs with perihelion in December (Pmax: 2222 ka, 2265 ka, and 2380 ka), red dots are runs with perihelion in June (Pmin: 2230 ka, 2346 ka, 2369 ka, and 2395 ka) -illustrating that eccentricity has a much larger effect on seasonality than precession at a given eccentricity.

	MDI (Log 10 (pg))	0.4 0.5 0.6 0.7			0.30 0.35 0.40	Seasonal contrast NPP (gC m -2 yr -1 )
		0.3
		0	0.02	0.04	0.06
			Eccentricity

  :The age model for Site MD05-293042 , in the Coral Sea, is based on ten AMS 14 C measurements of G. ruber sensu stricto (ss), and on correlation of the MD05-2930 δ 18 O G. ruber ss record with the LR04 benthic foraminiferal δ 18 O stack18 .IODP Site U1443[START_REF] Clemens | & the Expedition 353 Scientists Indian monsoon rainfall[END_REF] : The age model for Site U1443, on the Ninetyeast Ridge in the southernmost Bay of Bengal, is based on correlation of physical properties data on the primary shipboard splice to the LR04 benthic foraminiferal δ 18 O stack 18 . Physical properties were measured onboard the RV Joides Resolution on archive halves from all holes at Site U1443, and correlation and splicing were based on magnetic susceptibility, natural gamma ray, and reflectance spectroscopy b* data[START_REF] Clemens | Proceedings of the International Ocean Discovery Program; Indian monsoon rainfall; Expedition 353 of the[END_REF] .IODP Site U1448[START_REF] Clemens | & the Expedition 353 Scientists Indian monsoon rainfall[END_REF] : The age model for Site U1448, in the Andaman Sea, is based on correlation of the Site U1448 Ti/Ca record, derived by X-Ray Fluorescence (XRF) core scanning, to the XRF Ti/Ca record of Indian National Gas Hydrate Program (NGHP) Site 17 -a very close site with age control based on a benthic foraminiferal oxygen isotope record tuned to the LR04 benthic foraminiferal δ 18 O stack18,45,46 .IODP Site U1485 47 : The age model of the northern Papua Margin Site U1485 is based on correlation between the benthic foraminiferal δ 18 O record of this site and the LR04 benthic foraminiferal δ 18 O stack18 .IODP Site U1486 47 : The age model of the upper 100 m (the last 1.6 Ma) of Site U1486, from the southern bank of Manus Island, is based on a correlation between the benthic foraminiferal δ 18 O record and the LR04 benthic foraminiferal δ 18 O stack 18 . Below 100 m, the age model is based on the shipboard integrated magneto-biostratigraphy and correlation of the U1486 L* record (that shows strong G/IG variability) with the LR04 benthic foraminiferal δ 18 O stack 18 .

	IODP Site U1446 43 : The age model for Mahanadi Basin Site U1446 is based on correlation of a high-resolution benthic foraminiferal oxygen
	isotope stratigraphy to the LR04 benthic foraminiferal δ 18 O stack 18 (S.C. Clemens, unpublished data, pers. comm.)
	IODP Site U1483

47

: The age model for Scott Plateau Site U1483 is based on the shipboard integrated Magneto-biostratigraphy, and on correlation of the U1483 L* record, which presents strong G/IG variability, to the LR04 benthic foraminiferal δ 18 O stack 18 .

Table 1 :

 1 also simulate enhancement of the surface ocean temperature cycle at high eccentricity in the Eastern Equatorial Pacific, with higher amplitude that in the Western Equatorial Pacific. Characteristics of the nine marine records used in this study. Biome codes 20 : INDE = Eastern India Coastal province, MONS = Indian Monsoon Gyres Province, WARM = Western Pacific Warm Pool province, SUND = Sunda-Arafura Shelves Province, AUSE = East Australian Coastal Province, ISSG = Indian South Subtropical Gyre Province. These cores were retrieved during four expeditions: IMAGES 3[START_REF] Chen | Exploring quaternary variability of the east Asia monsoon, Kuroshio Current, and western Pacific warm pool systems: High-resolution investigations of paleoceanography from the IMAGES III (MD106) IPHIS cruise[END_REF] and 13[START_REF] Beaufort | Shipboard Scientific Party[END_REF] , and International Ocean Discovery Program (IODP) Expeditions 353[START_REF] Clemens | & the Expedition 353 Scientists Indian monsoon rainfall[END_REF] and 363 47 , which took place in the Western Pacific and Indian Tropical Oceans (Figure1).

	Extended Data Figures and Tables
	Site	Expedition	Longitude	Latitude	Water	Biome Record	Mean	Number of	Number of
						Depth		length		Resolution	samples	Measurements
						(m)		(kyr)		(kyr)
	U1446	IODP 353	85°441		19°050 N	1440	INDE	468		0.5	952	222 119
	U1448	IODP 353	93°00.00	10°38.03 N 1096	MONS 716		1.1	624	617 138
	U1443	IODP 353	90°21.71	5°23.01 N	2940	MONS 2818		1.6	1717	1 337 793
	MD97-2140	IMAGES 3	141°45.49	2°02.59 N	2547	WAR	1758		1.4	1241	540 523
							M
	U1486	IODP 363	144°36.08	2°22.34 S	1332	WAR	2661		2.2	1195	1 004 860
							M
	MD05-2920	IMAGES 13 144°32		2°51.48 S	1843	WAR	387		1.1	353	436 984
							M
	U1485	IODP 363	142°47.59	3°06.16 S	1145	SUND 535		0.9	587	580 612
	MD05-2930	IMAGES 13 146°15.73	10°25.30 S 1490	AUSE 793		2.1	372	620 736
	U1483	IODP 363	121°5.24	13°05.24 S 1733	ISSG	2332		2.3	994	1 777 406
										Average = 1.6 Sum =8035	Sum=7 138 171
	Extended Data Configuration (kyr)	Eccentricity Precession	Longitude of perihelion (°)	Obliquity (°)	Mean YCR at Equator (W.m -2 )
		C2222 (EmaxPmax) 0.054	0.054		87	23.73	74
		C2230 (EmaxPmin)	0.053	-0.037	315	23.84	83
		C2265		0.026	0.024		115	24.28	59
		C2346		0.034	-0.033	287	24.27	62
		C2369		0.016	-0.015	256	22.26	52
		C2380 (EminPmax)	0.006	0.006		77	23.74	37
		C2395 (EminPmin)	0.005	-0.003	310	23.42	39
	Extended Data

Table 2 :

 2 Summary of orbital parameters 19 used for each simulation and Mean Yearly Contrast of Radiation at Equator (Wm -2 ) derived from IPSL-CM5A2 (see Methods).

  Schematic illustrating how MDI is calculated (also see Methods), and how it represents morphological diversity, using data from Site U1485 as an example. In each sample, two size classes are created: coccoliths shorter and longer than 3 µm (blue and red histograms). The average log(mass) is calculated in both classes (blue and red dots on the histograms). MDI is the difference between these two

	a -Cross BT : Ecc. vs Size stack
		0.8		405 kyr	126 kyr	93 kyr
	Coherency	0.4 0.6					CL 90% CL 80%
		0.2
	Phase (rad)	-4 -2 0 2 4	0.92 Rad = 60 kyr Bandwidth	0.93 Rad = 18 kyr
		0		0.005		0.010
					Frequency
	b -MDI U1443		Period (kyr)			c -MDI U1486	Period (kyr)			d -MDI U1483	Period (kyr)
		405 126 96	41		23	19	405 126 96	41	23		19	405 126 96	41	23	19
	(Ka)							(Ka)					(Ka)
	Age							Age					Age
	MDI								MDI				MDI
	e -MDI MD97-2140 405 126 96 Extended Data Figure 2:	Period (kyr) 41	23	19	f -MDI MD05-2930 405 126 96	Period (kyr) 41	23	19	405 126 96 g -MDI U1448	41 Period (kyr)	23	19
	Age (Ka)							Age (Ka)				Age (Ka)
	MDI								MDI				MDI
	h -MDI U1485		Period (kyr)			i -MDI U1446 405 126 96	Period (kyr) 41	23	19	j -MDI MD05-2920
		405 126 96		41	23	19
													Period (kyr)
													405 126	96	41	23	19
	(Ka)								(Ka)
	Age								Age
	MDI								MDI				Frequency (cycle /kyr)

averages. Thus, MDI can differ substantially from average coccolith mass of the whole population (black dots on mass distribution plots), and two samples can have a similar MDI but very different average mass and mass distributions (e.g., examples at 180 ka and 190 ka), or very different MDI but similar averages (e.g., 3 ka vs 190 ka). In this way, MDI isolates morphology variability that results from evolutionary contractions and dilations in the number of different morphotypes present within the population.

Beaufort, L. et al. Insolation cycles as a major control of equatorial Indian ocean primary production. Science 278, 1451-1454, doi: 10.1126/science.278.5342.1451 (1997).

Thierstein, H. R., Geitzenauer, K. R., Molfino B. & Shackleton, N. J. Global synchroneity of late Quaternary coccolith datum levels: validation by oxygen isotopes. Geology 5, 400-404, doi: 10.1130/0091-7613(1977)5<400:GSOLQC>2.0.CO;2 (1977).

Raffi, I. et al. A review of calcareous nannofossil astrobiochronology encompassing the past 25 million years. Quaternary Sci. Rev. 25, 3113-3137, doi:https://doi.org/10.1016/j.quascirev.2006.07.007 (2006).

Longhurst, A. Ecological Geography of the Sea. Vol. 1 (Academic Press, 1998).

Acknowledgments

This paper is a contribution of the Climate research group at CEREGE. This research uses samples provided by the IODP, which is sponsored by the U.S. National Science Foundation and participating countries under management of Joint Oceanographic Institutions, Inc. We thank the shipboard participants of IODP Expeditions 353 and 363 and IMAGES Expeditions 3 and 13. The authors acknowledge French ANR projects CALHIS (LB), iMonsoon (CTB), AMOR (YD), and INSU project CALVE (CTB) which provided funding for this work. IODP-France provided post-cruise funding to LB and CTB. A. Fruy and S. Sergi are thanked for sample preparation assistance. We thank the CEA/CCRT for providing access to the HPC resources of TGCC under the allocation 2019-A0070102212 made by GENCI.

Author Contributions

LB designed the study. LB, YG, NB and MT developed automated artificial intelligence methods. LB, CTB, JCM, PC, EG, and SB prepared samples and/or generated data. ACS designed and ran the model simulations, in collaboration with YD. LB and CTB analysed the morphometric data. LB, CTB, ACS, BSM, YD, and YR discussed interpretations. LB, CTB and ACS wrote the manuscript with important contributions from BSM, YD, and YR.

Additional Information

Correspondence and requests for materials should be addressed to Luc Beaufort (beaufort@cerege.fr). Reprints and permissions information is available at www.nature.com/reprints.

Competing interest declaration

The authors declare no competing interests.

Data and Code availability

All morphological data are available at www.microautomate.cerege.fr/data/evolseason. 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 -OASIS3-MCT 2.0_branch (rev 4775 IPSL server)

The login/password combination requested at first use to download the ORCHIDEE 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.


Thomas Leduc
email: thomas.leduc@crenau.archi.fr

Vincent Tourre

Myriam Servières

The space-time cube as an effective way of representing and analysing the streetscape along a pedestrian route in an urban environment

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction and state of the art 1.Urban space analysis and visibility issues

Because a strong trend at work today is to develop and design cities that allow people to adapt and mitigate phenomena related to major issues, such as climate change, several new criteria have emerged for the design of sustainable and healthy cities. These new criteria have to do with walkability, role of the streetscape in encouraging the development of pedestrian mobility, physical activity and well-being. As [START_REF] Yin | Street level urban design qualities for walkability: Combining 2D and 3D GIS measures[END_REF] points out, indicators related to walkability are not only the D variables (design, density, diversity, destination accessibility, or distance to transit) and must necessarily include "the microscale streetscape features that affect the pedestrian experience."

This attention to the landscape and its influence on pedestrian behaviour is part of a long line of work on the interactions between individuals and their built environment in terms of environmental psychology. It led to the identification of about fifty urban qualities likely to influence the behaviour of pedestrians moving in the street space [START_REF] Ewing | Identifying and Measuring Urban Design Qualities Related to Walkability[END_REF].

The literature on visibility studies extends this questioning and provides definitions and operating procedures; tools that make it possible to formalize and to measure the essence of places in terms of some of the involved urban design qualities. In the state of the art, which aims to "explore the spatial properties of the visual structure inherent in space" [START_REF] Llobera | Extending GIS-based visual analysis: the concept of visualscapes[END_REF], we can highlight the isovist (Figure 1a) and isovist field (Figure 1b) as pioneering methodological tools proposed by [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF]. The isovist is indeed a tool that easily answers the questions raised by [START_REF] Batty | Exploring isovist fields: space and shape in architectural and urban morphology[END_REF]: "How far can we see, and how much can we see?". All key questions according to the author "in developing good design".

Even if other tools exist, such as the sky opening studied by [START_REF] Teller | A spherical metric for the field-oriented analysis of complex urban open spaces[END_REF], to analyse and even design urban space, we will focus in the current article only on isovists for their ability to synthesize, at any generation point, the immediately surrounding open space. The darker the colour, the more the corresponding portion of urban space is seen along the path.

Visibility representation issues along an urban route

According to [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF]'s definition, in an environment composed of opaque facades, the isovist associated with a given vantage point corresponds to all the points seen from the corresponding position in all directions of space. In the case of a 2D planar environment, such as the one presented in [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF]'s pioneering article, the isovist is a bounded polygon representing a horizontal slice of the pedestrian's immediate surrounding space. In the absence of any opaque surface, it may be useful to limit the scope of the gaze with an artificial horizon. This location-specific visual shape holds the same characteristics as a "visual pattern" as defined by [START_REF] Steiniger | Enabling pattern-aware automated map generalization[END_REF] and therefore information extracted from the isovist provides a synthetic view of the studied terrain.

With reference to the work on Gibson's motion perspective, [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF] proposes "a cumulative understanding of the form of the environment […] by perceiving variants and invariants in the transformation […] caused by the observer's movement." Indeed, because "the view of and visual exposure to the environment will gradually and sometimes suddenly change with [the] position", the corresponding set of isovists and their shapes and sizes, but also the nature of their boundaries, vary along the way so as to create some sort of footprints' field specific to a given environment. To move beyond the study of single isovist and to provide an analytical tool to study a set of positions, he therefore proposes the isovist field.

However, this cumulative visual experience exposes the designer to a difficulty in terms of representation, as can be seen in Figure 1b. The juxtaposition of polygons in the same horizontal plane causes overlaps. To overcome this difficulty and exhibit the structure of the set of surrounding visual spaces, [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF] proposes two strategies. The first consists in moving a panoptic light source along the path in an analogue mock-up. All the illuminated points of the mock-up correspond to the 3D isovist associated with the position of the light source. The second, described as more analytical by its author, consists in building a Minkowski model. This model is a sequential stack of 2D isovists where the Z-axis embeds the time dimension. Finally, the author proposes to apply a dimensional reduction to this isovist field and to study scalar fields (the isovist is abstracted at any point by a measure value such as its area). This simplified representation model makes it possible to derive iso-potential lines in a traditional way.

The space-time cube as a framework for temporal data visualization

A pedestrian's route through the city is a complex object that has a temporal dimension in addition to traditional spatial dimensions. As such, it is time-geography insofar as time and space are fundamental dimensions for the analysis of its dynamics. [START_REF] Hägerstrand | What about people in Regional Science?[END_REF] proposed innovative visual objects such as the space-time cube to illustrate in a graphic way this type of phenomenon. More recently, [START_REF] Bach | A Descriptive Framework for Temporal Data Visualizations Based on Generalized Space-Time Cubes[END_REF] has developed a conceptual framework for the visualization of temporal data based on generalized space-time cubes. According to [START_REF] Kraak | The Space-Time Cube Revisited from a Geovisualization Perspective[END_REF], "maps are used to stimulate visual thinking about geospatial patterns, relationships, and trends." The question is therefore whether the space-time cube is a relevant geovisual solution for the analysis of complex data such as an isovist field.

In this article, after a presentation of the conceptual framework developed by [START_REF] Bach | A Descriptive Framework for Temporal Data Visualizations Based on Generalized Space-Time Cubes[END_REF] and a review of some of the elementary operations resulting from their study, we will show that the representation techniques already implemented in the field of visibility analysis using isovists, fit perfectly into this formalism. With the help of this formalism, we will then propose innovative visualizations before opening up to extended visualization challenges.

Materials and method

Conceptual representation of spatial data based on the space-time cube

According to [START_REF] Bach | A Descriptive Framework for Temporal Data Visualizations Based on Generalized Space-Time Cubes[END_REF], "a space-time object is a geometrical object within a space-time coordinate system". The latter is a three-dimensional Euclidean space that adds a temporal dimension to a two-dimensional geographical space. In such a space, it is very common to represent the trajectory of a moving point object in motion on a plane map background by a static 3D curve. However, in the specific case of the representation of an isovist field, it is a matter of displaying not a point trajectory, but a collection of polygons as stated before. The issue lies in adapting a process for visualizing the evolution of a dimensionless object over time, to the visualization of the evolution and superimposition of a 2D shape over time. Regarding the issue of the representation device, the space-time cube is by definition a 3D object and its 2D visualization is part of the 3D rendering problems (with occlusion, depth ambiguity, perspective distortion, colour distortion and navigation difficulties [START_REF] Bach | A Descriptive Framework for Temporal Data Visualizations Based on Generalized Space-Time Cubes[END_REF]). To bypass the edge effects of a projection of a space-time cube on a 2D support as computer screens, various techniques exist to reduce its complexity and dimensionality.

Review of some elementary operations related to the space-time cube

According to [START_REF] Bach | A Descriptive Framework for Temporal Data Visualizations Based on Generalized Space-Time Cubes[END_REF], "the ultimate goal of space-time cube operations is to transform a spacetime cube into a space-time object whose shape is compatible with the shape of the media employed to convey the information". With this constraint, they then develop a taxonomy of elementary operations that they group into four categories: extraction (selecting a subset of the space-time object), flattening (dimensional reduction by aggregation), geometry transformation, and content transformation.

In the extraction category, we distinguish the so-called time cutting operation. This operation consists in extracting a 2D subset of the space-time object that corresponds to a single time cut (slice without thickness). In the flattening category, we can mention the time flattening which is a dimensional reduction operation by aggregating time. In the category of geometric transformations, we can mention the time scaling operation, which aims to resize the space-time object along the temporal axis (an acceleration of the phenomenon under study is similar to a contraction and a slowing down to a dilation). Finally, in the category of content transformations, we can mention the time colouring operation, which aims to alter the colours of each time slice according to time.

Focus on compound operations formalism

The inherent complexity of certain space-time phenomena justifies the application of several elementary operations, in sequence, to the data set to be displayed. To formalize these complex processes, [START_REF] Bach | A Descriptive Framework for Temporal Data Visualizations Based on Generalized Space-Time Cubes[END_REF] introduced the notion of compound operation. Like the usual composition of functions in mathematics, a compound operation is a mechanism to combine elementary operations to build composite ones. To develop this formalism, the authors introduced in particular the operator + which corresponds to the composition and the operator * which refers to a repetition (operation applied several times).

Discrete time flattening, which is an operation that aims to reduce the complexity of the space-time cube by aggregating all data independently of the time-related dimension, is then written as follows: (time cutting)* + time flattening. It consists of starting with a first step of sampling the space-time cube as a function of time. Then to apply an aggregation operation on the resulting layers while preserving the spatial dimensions.

Case Study

Presentation of the study site

In order to test the few techniques for representing the above-mentioned complex geographical objects, we have identified an urban fabric that is emblematic of the city of Nantes (France). This 700 m by 500 m portion offers a diversity of urban morphologies, on a mainly pedestrian sector, flat ground of constant altimetry. To survey this urban fragment, we have chosen a 1.6 km pedestrian route (see Figure 1b) that winds its way from the medieval district in the east to a district of classical architecture dating back to the 18 th and 19 th centuries in the west. A research study conducted in 2010 on this same route [START_REF] Leduc | Towards a spatial semantics to analyze the visual dynamics of the pedestrian mobility in the urban fabric[END_REF] showed that its direction, if it is perfectly explained from a chronological point of view, is perfectible in terms of better visibility of heritage objects. For obvious reasons related to the constrained size of this article, the illustrations we present from now on are limited to a small portion of this path (from sampling point #30 to sampling point #40, see Figure 1b).

Review of techniques already implemented for 2D visualization of isovists

The isovist presented in Figure 1a is a flat representation of the space-time cube at a given time step. It results from a single elementary time cutting operation. The isovist field shown in Figure 1b is a flat representation of the initial space-time cube. It is the result of a compound operation namely discrete time flattening described in section 2.3. To be more precise, this compound operation embeds a content transformation operation (translucent filter) in order to produce an alpha blending aggregation. The resulting compound operation is written: (time cutting + time colouring)* + time flattening. The same isovist field could have been shown with several isovists side-by-side as in [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF], thus corresponding precisely to the discrete time juxtaposing compound operation. It is indeed a matter of extracting "multiple time slices and placing them side-by-side or on a grid" [START_REF] Bach | A Descriptive Framework for Temporal Data Visualizations Based on Generalized Space-Time Cubes[END_REF].

The Minkowski model mentioned in [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF]'s article is a 3D rendering resulting from a compound operation applied to the space-time cube of isovist field: (time cutting + space shifting)* + 3D Rendering. In case this sequential stack of 2D isovists is coloured, as in Figure 2, an additional content transformation is applied and the corresponding compound operation is equal to (time cutting + space shifting + time colouring)* + time scaling + 3D Rendering. The time scaling operation that precedes the final 3D rendering step aims to facilitate the overall visualization of the shape by contracting it (by a factor of 2).

Fig. 2. Transposition of the "Minkowski model of the isovists along [the]

given path" proposed by [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF]. Sixteen isovists are stacked on the time axis, which is represented by the blue vertical axis. All the corresponding generation points are aligned on this axis. The more yellow the polygon is, the further away from the starting date.

The "Analog production of isovists along a path by point-source illumination of a model" stated in [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF] is not a static visualization but a dynamic operation that consists of moving the cutting position over time. This is the compound operation known as animated time cutting according to the above-mentioned formalism.

Innovative 3D renderings

The stack of isovists presented in Figure 2 has the advantage of providing a synthetic view of the evolution of the pedestrian's isovists during his walk in the city. Its formal qualities make it possible to highlight the bifurcations and configuration changes in the pedestrian's surroundings over time. However, the space shifting operation that systematically accompanies time cutting results in a loss of all anchor points because it translates all isovists into the same vertical.

To counteract this loss of anchoring and avoid getting lost on the background map, we propose to remove the space-shifting step in the aforementioned compound operation. The space-time object presented in Figure 3 is a 3D rendering resulting from the following compound operation applied to the generalized space-time cube: (time cutting + time colouring)* + time scaling + 3D Rendering. Fig. 3. By removing the space shifting operations in the Minkowski model from Figure 2, the isovists are relocated to match the path in spatial dimensions.

The shape presented in Figure 3 facilitates, through its stepped approach, a joint spatial and temporal interpretation of the isovist field. However, even if its sampling step is a little finer than that of Figure 1b (12 m vs. 20 m), one must recognize that this discretization remains coarse and that it may mask possible major visual events. Figure 4 combines a response in terms of refinement of discretization, the sampling step is reduced to 2 m, with a thickness representation of the corresponding "visual volume". Indeed, instead of a stack of isovist slices, we represent a global volume by connecting the contour points from one time step to the next time step. This rendering also has the effect of guaranteeing a form of continuity of the entire shape, that could be associated to a single 3D Rendering operation, and thus allowing to show the raw shape of the isovist field. Fig. 4. Rendering in thickness of the space-time shape presented in Figure 3. In addition to refining the sampling step, the contour surfaces are drawn to shape a kind of "visual volume".

Discussion

Immersion

This article and the figures it embeds are intrinsically 2D objects. Being able to free ourselves from the 3D rendering step, by implementing stereographic rendering, would probably allow us a better understanding of an isovist field along a path (which is a typically 3D object). Indeed, by avoiding projection on a plane, we give free opportunity to immersion and navigation in the virtual space-time cube. We also allow to resize and adapt it to a human scale, as [START_REF] Bigot-Doll | Échelles diffuses : Prospectives de conception via matière et VR[END_REF] points out: "It can be seen that the scanned object […] is systematically enlarged by students within the virtual environment in order to reach a height at least equal to the size of the user. The model is now 'inclusive' and in turn operates a form of immersion that we will call expansive".

Increase dimensionality

As usual in the representation of the space-time cube, the background map in Figures 2,3, and 4 is two-dimensional. To facilitate the location and identification of built volumes, we could imagine replacing this background plan with a three-dimensional urban model. However, this facility would have the disadvantage of introducing new masks and thus increasing occlusions. It would also introduce confusion in the use of the third dimension, which would be both spatial and temporal. In order to analyse the interstitial urban space in a volumetric perspective, we could imagine representing not two-dimensional isovists in their evolution over time but three-dimensional isovists. This dimensional increase of the item to be represented will consequently substantially increase the occlusion effects and depth ambiguities.

About adaptive sampling

The sampling distances used in this study are arbitrary and systematic. However, there are portions of the route that present more or less visual events. It is therefore entirely legitimate, if not desirable, to consider adapting the sampling step to the regularity or irregularity of the course at each point. The partitioning of the urban space of [START_REF] Peponis | On the description of shape and spatial configuration inside buildings: convex partitions and their local properties[END_REF], as well as its extension by [START_REF] Leduc | Motion Perspectives Integration in the Qualification of the Urban Spaces: Towards a 2D-and a 3D-Enrichment of the S-Partition Method[END_REF], who proposes to qualify the space according to the importance of changes in visual surroundings, could be very useful in this regard.

Conclusion

In this article, we have stated that attention to the visual landscape surrounding pedestrians as they walk around the city is a useful operating concept in urban planning processes. To this end, [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF] proposed the isovist field as a solution allowing a dynamic analysis (spatial and temporal) of the space covered. The representation of a complex polygonal object, such as a 2D isovist, in its evolution over time along a pedestrian path, is however a hard research issue. From the generalized space-time cube proposed by [START_REF] Bach | A Descriptive Framework for Temporal Data Visualizations Based on Generalized Space-Time Cubes[END_REF], which is a conceptual framework adapted to the analysis and representation of the isovist field, we have resituated all the representation solutions that have already been proposed. In this article, we have also identified the problems raised by each of these solutions and proposed innovative evolutions. We now have to validate the usability of our proposals through a set of controlled experiments.

Fig. 1 .

 1 Fig. 1. (a) On the horizontal plane of the gaze, the isovist is the set of points seen from a given position (black disc); (b) The 1.6 km long pedestrian route is sampled every 20 m. An isovist is then calculated at each sampling point. The corresponding isovist field is drawn in a flat representation. The darker the colour, the more the corresponding portion of urban space is seen along the path.

SHS Web of Conferences 64, 03005 (2019) https://doi.org/10.1051/shsconf/20196403005

SHS Web of Conferences 64, 03005 (2019) https://doi.org/10.1051/shsconf/20196403005