A 15-million-year-long record of phenotypic evolution in the heavily calcified coccolithophore Helicosphaera and its biogeochemical implications

The biogeochemical impact of coccolithophores is defined not only by their overall abundance in the oceans but also by wide ranges in physiological traits such as cell size, degree of calcification and carbon production rates between different species. Species' sensitivity to environmental forcing has been suggested to relate to their cellular PIC : POC (particulate inorganic carbon : particulate organic carbon) ratio and other physiological constraints. Understanding both the short-term and longer-term adaptive strategies of different coccolithophore lineages, and how these in turn shape the biogeochemical role of the group, is therefore crucial for modeling the ongoing changes in the global carbon cycle. Here we present data on the phenotypic evolution of a large and heavily calcified genus Helicosphaera (order Zygodiscales) over the past 15 million years (Myr), at two deep-sea drill sites in the tropical Indian Ocean and temperate South Atlantic. The modern species Helicosphaera carteri, which displays ecophysiological adaptations in modern strains, was used to benchmark the use of its coccolith morphology as a physiological proxy in the fossil record.

Our results show that, on the single-genotype level, coccolith morphology has no correlation with growth rates, cell size or PIC and POC production rates in H. carteri. However, significant correlations of coccolith morphometric parameters with cell size and physiological rates do emerge once multiple genotypes or closely related lineages are pooled together. Using this insight, we interpret the phenotypic evolution in Helicosphaera as a global, resource-limitation-driven selection for smaller cells, which appears to be a common adaptive trait among different coccolithophore lineages, from the warm and high-CO₂ world of the middle Miocene to the cooler and low-CO₂ conditions of the Pleistocene. However, despite a significant decrease in mean coccolith size and cell size, Helicosphaera kept a relatively stable PIC : POC ratio (as inferred from the coccolith aspect ratio) and thus highly conservative biogeochemical output on the cellular level. We argue that this supports its status as an obligate calcifier, like other large and heavily calcified genera such as Calcidiscus and Coccolithus, and that other adaptive strategies, beyond size adaptation, must support the persistent, albeit less abundant, occurrence of these taxa. This is in stark contrast with the ancestral lineage of Emiliania and Gephyrocapsa, which not only decreased in mean size but also displayed much higher phenotypic plasticity in their degree of calcification while becoming globally more dominant in plankton communities.

1 Introduction

Coccolithophores (calcifying haptophyte algae) are a globally abundant group of marine phytoplankton and an important component of the biogeochemical carbon cycle. Their calcification (Monteiro et al., 2016), primary production (Rousseaux and Gregg, 2013) and deep-sea burial of calcite scales (coccoliths; Baumann et al., 2004; Ziveri et al., 2007) have been crucial for the evolution of past and modern ocean chemistry, and coccoliths are still alleviating the negative effects of rising atmospheric and oceanic CO₂ levels by the ballasting of organic carbon (Ridgwell and Zeebe, 2005). The cellular biogeochemical output of coccolithophores, which is commonly summarized as a cellular balance of particulate inorganic carbon (PIC) and particulate organic carbon (POC) production rates (i.e., PIC : POC ratio), is highly susceptible to environmental forcing such as high-temperature stress (Gerecht et al., 2018, 2014; Rosas-Navarro et al., 2016), nutrient limitation (Bolton and Stoll, 2013; Gerecht et al., 2015; Müller et al., 2017) and ocean acidification (Ridgwell et al., 2009; Riebesell et al., 2000). Understanding how these marine algae adapt to environmental changes and how adaptive strategies shape their biogeochemical impact is thus essential for modeling the carbon cycle dynamics in the projected warmer, high-CO₂ oceans of the near future (Bopp et al., 2013; Doney et al., 2004; Feely et al., 2004).

Deep-sea sediments hold a ∼200 Myr long fossil record of coccoliths (Bown et al., 2004). Some morphometric features of coccoliths, such as length and thickness, correlate with biogeochemically relevant traits: cell size (Henderiks, 2008), growth rates (Gibbs et al., 2013) and the production rates of PIC and POC (Bolton et al., 2016; McClelland et al., 2016). Combined with morphospecies diversity and community composition that is readily inferred from the fossil record (Bown et al., 2004; Knappertsbusch, 2000; Suchéras-Marx and Henderiks, 2014), coccolith morphometry could therefore help elucidate evolutionary processes from the cellular to the community level. In this context, the transition from the warm and high-CO₂ world of the middle Miocene (∼15 Ma) to the colder, low-CO₂ conditions of the Pleistocene (Herbert et al., 2016; Sosdian et al., 2018; Super et al., 2018; Zachos et al., 2001; Zhang et al., 2013) provides an optimal case study for investigating the long-term evolutionary adaptation of modern coccolithophore lineages to relevant climate forcing and understanding the biogeochemical implications of their evolutionary adaptation. This was a period of large-scale shifts in coccolithophore community structure, abundance and biogeochemical impact (Hannisdal et al., 2012; Si and Rosenthal, 2019), as also evidenced by a macroevolutionary decrease in cell size (Aubry, 2009; Bolton et al., 2016; Imai et al., 2015; Knappertsbusch, 2000; Suchéras-Marx and Henderiks, 2014; Young, 1990) and degree of calcification (Bolton et al., 2016).

To date, studies of coccolithophore evolutionary adaptation have mostly focused on members of the Noelaerhabdaceae family, which dominated coccolithophore communities during the middle Miocene to Pleistocene (e.g., Henderiks and Pagani, 2007; Suchéras-Marx and Henderiks, 2014). While their most important modern representative Emiliania huxleyi is highly abundant in today's oceans, this lineage represents only a small fraction of the morphological, ecological and physiological diversity of coccolithophores (Young et al., 2005). The Noelaerhabdaceae also have a unique evolutionary history with an early Jurassic divergence (∼195 Ma) from the other calcifying haptophytes (Medlin et al., 2008). Clearly, there is a need to investigate evolutionary adaptation across the range of diversity and include other prominent coccolithophore lineages to determine whether the evolutionary patterns observed in the Noelaerhabdaceae are universal for the group as a whole.

This study focuses on the phenotypic evolution of the heavily calcified coccolithophore genus Helicosphaera (order Zygodiscales; Fig. 1). This genus was a constant component of coccolithophore communities during the Cenozoic (Aubry, 1988) and is still an important contributor to oceanic carbonate fluxes today (Baumann et al., 2004; Menschel et al., 2016; Ziveri et al., 2007). The most common modern species within this genus, Helicosphaera carteri (Wallich) Kamptner, can be grown in culture (Sheward et al., 2017; Šupraha et al., 2015) and has a long fossil record spanning over 20 Myr (Aubry, 1988). Earlier experimental studies have found a range of strain-specific physiological traits in H. carteri, likely related to local ecophysiological adaptation (Šupraha et al., 2015). In addition, this species is considered highly susceptible to ocean acidification due to its heavily calcified coccoliths and high PIC : POC ratio (Gafar et al., 2019). Finally, both modern and fossil Helicosphaera morphospecies are well-defined and identifiable under a light microscope, which allows for long-term tracking of speciation events, community dynamics and phenotypic evolution at the morphospecies level. The main aims of this study are to (i) benchmark the use of coccolith morphology as a potential physiological proxy in Helicosphaera based on the coccolith morphology of physiologically distinct H. carteri strains in culture, (ii) explore the adaptive strategies of the genus and main evolutionary drivers over the past 15 Myr, and (iii) compare the evolutionary patterns of Helicosphaera with other major coccolithophore lineages and discuss the biogeochemical implications of their evolutionary adaptation.

Figure 1 Helicosphaera coccosphere and coccolith morphology. (a) Scanning electron micrograph of a Helicosphaera carteri coccosphere and a polarized light microscope image of an individual coccolith obtained from culture (Atlantic strain RCC 1323). (b) Light microscope images of coccoliths belonging to fossil Helicosphaera species analyzed in this study: H. carteri, H. granulata and H. sellii. Scale bar = 2 µm. (c) Schematic illustration of Helicosphaera coccolith morphology and calcite birefringence under cross-polarized light (after Young et al., 2004). Note the distinction between the birefringent blanket area (white) and the partially birefringent area of the coccolith wing or flange (grey).

2 Material and methods

2.1 Helicosphaera carteri cultures

Culture data, obtained from batch experiments presented in Šupraha et al. (2015), are herein used as a benchmark for phenotypic plasticity in modern H. carteri. Two strains were isolated from contrasting environments, the eutrophic South Atlantic Ocean (strain RCC1323) and oligotrophic western Mediterranean Sea (strain RCC1334; Fig. 2a). These strains exhibited different physiological traits in terms of, for example, cell size, growth and calcification rates (PIC production rates). Strains were acclimated to experimental conditions for 10 generations prior to inoculation. Experiments were run in triplicate at a temperature of 17 ^∘C and an irradiance of ∼160 µmol m⁻² s⁻¹ under a 14 : 10 h light : dark cycle. Further details of the experimental setup are given in Šupraha et al. (2015). All phenotypic and physiological data used in the present study were obtained from the exponentially growing, nutrient-replete control experiments only.

Figure 2 (a) Locations of DSDP Site 525 and ODP Site 707 (stars) and other coring sites (open circles) under discussion, as well as isolation sites of RCC cultures (squares) analyzed in Šupraha et al. (2015). Map as Mollweide projection with an overlay of annual average Chlorophyll a concentrations (data obtained from the Ocean Color database, NASA, 2018, http://oceancolor.gsfc.nasa.gov/, last access: 2 August 2019) was created using the Ocean Data View 4 software (Schlitzer, 2018). (b) Long-term trends in sea surface temperature (SST) at selected temperate and tropical sites (data from Herbert et al., 2016).

2.2 Deep-sea sediment records

Fossil time series data were collected from pelagic foraminiferal nannofossil oozes and chalks recovered at two deep-sea drill sites (Fig. 2a). Deep Sea Drilling Project (DSDP) Leg 74, Site 525, is located in the southeastern Atlantic Ocean, at the Walvis Ridge (29^∘04.24^′ S, 02^∘59.12^′ E; water depth 2467 m; Moore et al., 1984). Ocean Drilling Program (ODP) Leg 115, Site 707, is located in the tropical western Indian Ocean, at the Mascarene Plateau (07^∘32.72^′ S, 59^∘01.01^′ E; water depth 1552 m; Backman et al., 1988). The sites were selected to represent contrasting productivity levels and climatic conditions. Modern-day primary production estimates are on average 50–100 g C m⁻² yr⁻¹ higher in the South Atlantic than in the equatorial Indian Ocean (Antoine et al., 1996; Beaufort et al., 1997; Fischer et al., 2000). Consistently higher carbonate mass accumulation rates at Site 525 compared to Site 707 (Suchéras-Marx and Henderiks, 2014) confirm consistent site-to-site carbonate production contrasts during the geological past. Regional differences in mean sea surface temperature extend back at least 10 million years (Fig. 2b). Age models are based on calcareous nannoplankton biostratigraphy (this study; Backman et al., 1988; Moore et al., 1984; Suchéras-Marx and Henderiks, 2014), updated to the most recent geological timescale (Gradstein et al., 2012; Table S1 in the Supplement). Linear sedimentation rates were estimated between age–depth tie points.

2.3 Coccolith morphometry

Helicosphaera carteri coccoliths harvested from triplicate cultures were filtered on cellulose nitrate filters (0.8 µm, Whatman) after dispersing the coccoliths with a Triton–NaOCl treatment (Paasche et al., 1996). Filters were dried at 60 ^∘C and mounted with Canada balsam (Merck, USA) on microscope slides, rendering filters transparent under polarized light. Microscope slides for fossil coccolith morphometry were prepared from fine fraction (<38 µm) sediment using the spraying technique (Henderiks and Törner, 2006) for a total of 36 (Site 525) and 38 (Site 707) stratigraphic levels. Subsamples of freeze-dried bulk sediments were wet-sieved with buffered distilled water into three separate size fractions (<38, 38–63 and >63 µm). The fine fraction, dominated by coccoliths, was collected in suspension, concentrated on cellulose acetate membrane filters using a Millipore vacuum pump and oven-dried at 50 ^∘C.

In total, 331 culture-harvested and 4115 fossil Helicosphaera specimens were analyzed using a Leica DM6000 B microscope equipped with polarizing filters placed at different angles (Beaufort et al., 2014). Besides size measurements, mean grey level (MGL) was used as a proxy for a change in coccolith thickness based on the birefringence of calcite (Beaufort, 2005). Standard imaging and light settings were ensured during the course of image collection (intensity = 200; field = 10; aperture = 10). Image collection was performed within a short time frame in order to minimize the effect of light bulb aging. For each sample, at least 50 Helicosphaera coccoliths were imaged at 1000× magnification with a SPOT Flex color digital camera (Diagnostic Instruments Inc.; micrometer-scale calibrated pixel resolution 0.061 µm (1.2× C-mount) or 0.074 µm (1.0× C-mount)) under three different orientations of the polarizer (0, 35 and 45^∘). Coccoliths were randomly selected and identified to the morphospecies level. Images taken with the three orientations of the polarizer were merged in a single 8 bit greyscale image, showing no extinction cross (Beaufort et al., 2014; Fig. 1). Images from each sample were analyzed as a batch using a custom-made macro in ImageJ software v1.47c. It determined the MGL of the selected area, measured in grey-level values of 0 (black) to 255 (white), and performed size measurements on individual coccoliths (area, perimeter, and length and width of the fitted ellipse). The calculation of coccolith thickness and volume was conducted following the protocol described in Beaufort et al. (2014):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{Thickness}\left(\mathrm{\mu }\mathrm{m}\right)=\mathrm{MGL}\cdot (\mathrm{1.55}/\mathrm{255})\left(\mathrm{\mu }\mathrm{m}\right),\text{(2)}& {\displaystyle }\mathrm{Volume}\left(\mathrm{\mu }{\mathrm{m}}^{\mathrm{3}}\right)=\mathrm{Thickness}\left(\mathrm{\mu }\mathrm{m}\right)\cdot \mathrm{Area}\left(\mathrm{\mu }{\mathrm{m}}^{\mathrm{2}}\right).\end{array}$$

It should be noted that all estimates of Helicosphaera coccolith volume are systematically underestimated but that this does not preclude accurate, relative comparisons between samples. Variations in area estimates are reliable within the constraints of image segmentation, but Helicosphaera coccoliths (also known as helicoliths) are not entirely birefringent under crossed or circular polarized light due to the crystallographic orientation of calcite elements; the proximal plate and the blanket that covers most of the central area are birefringent in analysis, but the flange that extends from the proximal plate and underlies the blanket is composed of so-called V-units that are only partially birefringent in the same orientation on a microscope slide (Fig. 1; Beaufort et al., 2014; Young et al., 2004). In addition, parts of the central area exceed 1.55 µm (or 255 (white) value; Beaufort et al., 2014), showing faint yellow birefringence colors that render grey values (<255) in 8 bit, thus lowering the mean grey-value estimate calculated from all pixels in one specimen. An examination of random specimens from both sites performed in ImageJ using the “Threshold color” function has shown that the yellow signal from the thickest parts of the central area usually accounts for up to 5 % of surface area in medium-sized species (e.g., H. carteri and H. sellii) and up to 10 %–15 % in the largest species H. granulata. Therefore, Helicosphaera volume estimates are not converted to calcite mass (multiplied by the density of calcite, 2.7 pg µm⁻³). Instead, in analogy with previous studies on Coccolithus pelagicus (Cubillos et al., 2012; Gerecht et al., 2014), we report on changes in the coccolith volume index (CVI; ln µm³) that scales with other coccolith size and shape parameters. Statistically significant changes in the CVI are meaningful indicators, if not a muted expression, of the true phenotypic changes in both modern and fossil Helicosphaera assemblages.

Finally, we compiled morphometric time series data for Calcidiscus (Knappertsbusch, 2000) and Reticulofenestra (Bolton et al., 2016; Hannisdal et al., 2012) from tropical and midlatitude locations to compare patterns of phenotypic variation and evolution across major coccolithophore clades. This comparison includes a total of 4907 Reticulofenestra coccoliths measured in 37 samples at Site 525. The latter database is limited to linear size measurements (length) only, because it was collected with a different polarized light microscopy and digital image analysis setup (Henderiks, 2008; Henderiks and Törner, 2006).

2.4 PIC : POC and geometric proxies for cell physiology

Cellular production rates of particulate inorganic carbon (PIC) and particulate organic carbon (POC) were determined for H. carteri in the exponentially growing control experiments as a product of growth rate and the corresponding cellular elemental quotas (Šupraha et al., 2015). The molar ratio PIC : POC is widely used to assess the relative balance between calcification and photosynthesis, which in turn informs us about the net CO₂ uptake or release by coccolithophores (review by Ridgwell et al., 2009, and references therein).

These physiological rates cannot be directly measured in the fossil record. Nevertheless, within a theoretical framework of the relationships between cell geometry, coccolith dimensions and coccolith coverage, cellular levels of organic and inorganic carbon (calcite) quotas can still be estimated with a range of approaches (Bolton et al., 2016; Gibbs et al., 2013, 2018; Henderiks, 2008; Jin et al., 2018; McClelland et al., 2016).

In this study, we use the approach described by McClelland et al. (2016), who correlated a coccolith shape index, the aspect ratio (AR), with PIC : POC ratios measured in modern cultured strains of Emiliania huxleyi and Gephyrocapsa oceanica, deriving the following power–law relationship:

$$\begin{array}{}\text{(3)}& \mathrm{PIC}:\mathrm{POC}=e^{\mathrm{3.5}\pm \mathrm{0.2}}{\mathrm{AR}}^{\mathrm{1.12}},\end{array}$$

where

$$\begin{array}{}\text{(4)}& \mathrm{AR}=\mathrm{Thickness}/\surd \left(\mathrm{Area}\right).\end{array}$$

Despite the fact that coccolith thickness in Helicosphaera is systematically underestimated, we took a similar approach to investigate relative changes in shape (AR) and cellular degree of calcification (PIC : POC) between mean phenotypes.

2.5 Coccolith abundance counts and paleofluxes

For the quantitative analysis of fossil coccolith assemblages, microscope slides were prepared directly from freeze-dried bulk sediments using the drop technique (Bordiga et al., 2015). The 74 samples were analyzed at 1000× magnification under an Olympus BX53 polarized-light microscope. A minimum of 300 coccoliths were counted and at least 20 fields of view (FOVs) were analyzed in each sample to determine the relative contribution of the major coccolithophore taxa. The absolute abundance of coccoliths (N g⁻¹) was calculated using the formula

$$\begin{array}{}\text{(5)}& A{A}_{\mathrm{coccolith}}=(N\times A)/(f\times n\times W),\end{array}$$

where N is the number of coccoliths counted, A is the area of the coverslip, f is the area of one FOV (0.0367 mm²), n is the number of FOVs counted and W is the weight of the sediment on the coverslip (g).

Fluxes of Helicosphaera spp. and other taxa (N m⁻² yr⁻¹) were calculated from the absolute abundances and the site-specific mass accumulation rates (MAR; data from Suchéras-Marx and Henderiks, 2014) with a reproducibility of ±15 %:

$$\begin{array}{}\text{(6)}& {\mathrm{Flux}}_{\mathrm{coccolith}}=A{A}_{\mathrm{coccolith}}\left(\mathrm{N}{\mathrm{g}}^{-\mathrm{1}}\right)\times \mathrm{MAR}\left(\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}\right).\end{array}$$

Past primary productivity (PP) estimates (g C m⁻² yr⁻¹) were calculated from the relative abundance of the lower-photic-zone species Florisphaera profunda (in percent, % Fp), taking the 95 % confidence intervals as upper and lower inputs (cf. Beaufort et al., 1997):

$$\begin{array}{}\text{(7)}& \mathrm{PP}=\mathrm{617}-[\mathrm{279}\times \log (\mathit{\%}\mathit{Fp}+\mathrm{3}\left)\right].\end{array}$$

The % Fp proxy is highly reliable for estimating PP at latitudes between 30^∘ N and 30^∘ S. At higher latitudes, and in upwelling regions, the relative abundance of F. profunda is also controlled by lower mean temperatures (Hernández-Almeida et al., 2019).

2.6 Statistical analysis

Raw biometry data (micrometer scale) were transformed to their natural logarithm (natural logarithmic micrometer scale) before calculating sample means and variances when investigating allometric relationships and comparing temporal trends within the time series.

One-tailed t tests were performed in R version 3.2.2. using the “t.test” command from the “Stats” package. Paired t tests were used to compare means between culture strains and between the fossil measurements. For the latter, site-to-site comparisons were made at the genus and (morpho)species level; p values < 0.05 were considered statistically significant. The 95% confidence intervals for relative abundance estimates were calculated in Past version 3.20 freeware (Hammer et al., 2001; Suchéras-Marx et al., 2019).

3 Results

3.1 Helicosphaera carteri allometry

Coccolith area and thickness show a positive correlation in H. carteri, although thickness scales up at a marginally slower rate than area (i.e., slightly negative allometry), which is highly similar for modern strains and fossil specimens (Fig. 3a, b). Mean coccolith dimensions in modern single genotypes fall within the lower range of the fossil measurements. Differences in mean coccolith dimensions between the South Atlantic and Mediterranean strains are not statistically significant or only borderline significant (paired t tests, p values = 0.029 to 0.039; Table S2). The range of shapes (AR) and dimensions (CVI) in cultured H. carteri largely overlap with its fossil counterparts (Fig. 3c, d). The Mediterranean strain (RCC1334), however, produced higher proportions of small and lightly calcified coccoliths (outliers from the distribution; Fig. 3d). Likewise, the raw fossil data ranges indicate that Indian Ocean Site 707 hosted higher proportions of small H. carteri coccoliths than South Atlantic Site 525 (Fig. S1, in the Supplement).

Figure 3 Coccolith allometry in modern strains and fossil specimens of H. carteri. (a–b) Natural log-transformed (ln) coccolith surface area (µm²) vs. thickness (µm) for all individual fossil specimens at sites 525 (blue; N=1156) and 707 (red; N=1399) and mean (±1 SD) values for South Atlantic (black square) and Mediterranean (open square) strains; (c–d) coccolith volume index (CVI; ln µm³) vs. aspect ratio (AR; unitless). Here, individual specimens are also shown for the modern strains.

3.2 Mean phenotypic evolution of the Helicosphaera genus

The fossil time series data, now considering all measurements of H.carteri and other fossil Helicosphaera morphospecies, show long-term trends in mean CVI and AR at both sites (Fig. 4a, b). Linear regression of sample means vs. age suggests a 16 % decrease in mean size and CVI (natural logarithmic scale) from the middle Miocene to late Pleistocene at South Atlantic Site 525 (r=0.849), which corresponds to a 46 % decrease in mean volume on the linear scale (µm³). The same trend is less pronounced at tropical Indian Ocean Site 707, with a 9 % decrease in mean size and CVI (27 % decrease in volume) over the same time interval (r=0.525). Mean shape (AR) estimates show an opposite trend with age at both sites, implying a 24 % (525; r=0.946) and 18 % (707; r=0.815) increase in aspect ratio over the past 15 million years. The AR is mainly driven by changes in size, since the decrease in thickness over time is less pronounced (Fig. S2), suggesting that the relative degree of calcification increased as coccolith dimensions (and thus mass per individual coccolith) decreased at both sites.

Figure 4 Long-term trends in phenotypic evolution of the Helicosphaera genus (including all morphospecies) and paleoenvironmental estimates at DSDP Site 525 (in blue) and ODP Site 707 (in red) over the past 15 Myr. (a) Helicosphaera mean coccolith volume index (CVI; ln µm³) and (b) coccolith aspect ratio (AR). N≈50 per sample. Error bars show 1 SD. (c) Primary productivity (PP) estimates (gC m⁻² yr⁻¹; 95 % confidence intervals) based on the relative abundance of Florisphaera profunda (% Fp); cf. Beaufort et al. (1997). Shading between 8 and 10 Ma indicates the interval of the first occurrence of this taxon. Bars to the left show the range in modern primary production estimates for the South Atlantic (blue) and the eastern Indian Ocean (red; Antoine et al., 1996; Beaufort et al., 1997; Fischer et al., 2000). (d) Long-term trends in placolith fluxes (or burial rates; N m⁻² yr⁻¹) of species dwelling in the upper photic zone only. Error bars show ±15 % reproducibility as determined from repeat analyses (Bordiga et al., 2015).

When comparing all fossil measurements over the past 15 Myr, mean CVI was significantly smaller at Site 707 than at Site 525 (paired t test; t=20.36, df = 3863.8, p value < 0.001). The Indian Ocean record also reveals larger variation in Miocene–Pliocene sample means compared to the more gradual trend observed in the South Atlantic, so site-to-site offsets in mean CVI varied over time (Fig. 4a). These offsets were largest throughout the middle Miocene, after which mean CVI reached more similar values that persisted during the late Miocene between ∼ 8.5 and 5.5 Ma. Site-to-site offsets increased again during the latest part of the Miocene and Pliocene, ∼ 5.5–3 Ma. From a low in the early Pleistocene (∼2.5 Ma) both sites display a near-synchronous increase in CVI that peaked during the middle Pleistocene (∼1.2 Ma). A steady decrease in CVI followed during the late Pleistocene at both sites; this trend appears steeper at Site 707, so mean CVI values diverged again between sites during the past ∼1 million years.

3.3 Phenotypic evolution at the morphospecies level

Both sites reveal distinct temporal patterns in Helicosphaera morphospecies composition (Figs. 5, S2). In total, nine morphospecies were observed: H. carteri, H. hyalina, H. wallichii, H. inversa, H. sellii, H. ampliaperta, H. granulata, H. euphrates and H. intermedia. Of these, three species (H. carteri, H. sellii and H. granulata; Fig. 1) were the most abundant, with common trends in relative abundances over the ∼15 Myr interval (Fig. 5a).

Figure 5 Phenotypic evolution of the most abundant Helicosphaera fossil morphospecies at DSDP Site 525 (in blue) and ODP Site 707 (in red) over the past 15 Myr. (a) Relative proportions (%) of H. carteri, H. sellii, H. granulata and other Helicosphaera morphospecies (closed sum). (b–d) Average CVI (ln µm³) for each of the three most prominent morphospecies (names indicated in plots). The number of measurements per fossil data point is between 2 and 50 depending on the relative abundance of the species in a sample.

Helicosphaera granulata was the most dominant species at the South Atlantic site during the middle Miocene (∼ 15–11 Ma) and occasionally within the same period (∼14.5 and ∼12 Ma) at Site 707. The rest of the interval was dominated by H. carteri, with higher contributions of H. sellii during the Pliocene and Pleistocene (∼ 4.5–1 Ma). Of the three studied species, only H. carteri was present continuously during the past 15 Myr (Figs. 5b, S1). Helicosphaera granulata was present until ∼ 2–3 Ma, whereas H. sellii had its first occurrence at ∼15 Ma at Site 525 and ∼12 Ma at Site 707 and its last occurrence at ∼1 Ma at both sites. Interestingly, extinctions of H. sellii and H. granulata were preceded by a decrease in their mean CVI at both sites (Fig. 5c, d). Similarly, H. carteri has experienced a distinct decline in the CVI over the past million years, and modern H. carteri strains have lower mean size and volume than most of their fossil sample counterparts (Figs. 3, 5b).

Consistently, the mean CVI of H. carteri, H. sellii and H. granulata was larger at Site 525 than at Site 707 (H. carteri t=10.83, df = 2427.2, p value < 0.001; H. sellii t=7.96, df = 324.91, p value < 0.001; H. granulata t=15.15, df = 890.57, p value < 0.001). Nevertheless, the long-term species-level patterns, as well as site-to-site offsets, were comparable among the three species (Fig. 5b–d), demonstrating that the genus-level coccolith volume trends and offsets between the sites (Fig. 4a) were caused not merely by relative contributions of differently sized species but also by phenotypic evolution at the species level.

3.4 Paleoecological records

Stark differences in nannofossil community composition and nannofossil fluxes reflect the contrasting environmental conditions between sites (Fig. S3). Site 525 had consistently higher nannofossil fluxes (average 6.49×10⁹ N m⁻² yr⁻¹) than Site 707 (average 2.98×10⁹ N m⁻² yr⁻¹) and was dominated by the opportunistic genera Reticulofenestra and Gephyrocapsa that bloom in the upper photic zone, while assemblages at Site 707 were dominated by the species Florisphaera profunda that is restricted to the deep-photic zone. Relative abundances of Helicosphaera spp. were relatively low (<10 % of assemblage), but it was consistently present at both sites with average fluxes of 7.57×10⁸ N m⁻² yr⁻¹ (Site 525) and 2.30×10⁸ N m⁻² yr⁻¹ (Site 707). Other heavily calcified taxa such as Coccolithus pelagicus s.l. (Coccolithaceae) had its highest relative abundance (5 %–15 % of assemblage) and fluxes (up to 1.92×10⁹ N m⁻² yr⁻¹) at Site 525 prior to ∼7 Ma.

Florisphaera profunda-based primary productivity estimates for the past 8–10 Myr (Fig. 4c; first occurrence of this taxon in Tortonian stage; http://www.mikrotax.org/Nannotax3, last access: 4 May 2020) and fluxes of upper-photic-zone taxa (total placolith fluxes in Fig. 4d) confirm that significant site-to-site differences in primary production were maintained throughout the past 15 million years. High % Fp (up to 75 %) at Site 707 indicates an oligotrophic system with low to moderate phytoplankton production in the upper photic zone (PP average 139±33 g C m⁻² yr⁻¹) compared to the more productive Site 525 where % Fp was consistently <20 % and, therefore, PP > 200 g C m⁻² yr⁻¹ (average 359±65 g C m⁻² yr⁻¹; Figs. 4c, S3). The highest nannofossil fluxes were recorded during the late Miocene and earliest part of the Pliocene (∼9 and 4 Ma), when lowered % Fp values also corroborate elevated PP levels at both sites. Despite the highly different community structures and productivity, a similar transition from larger Reticulofenestra species that dominated the middle Miocene community to smaller Gephyrocapsa spp. in the late Pliocene and Pleistocene is observed at both sites (Fig. S4).

4 Discussion

4.1 Are coccolith dimensions a proxy for cell physiology?

Mean coccolith dimensions (CVI, AR or any other parameters tested) in modern Helicosphaera carteri strains did not reflect the significant differences in growth rates, cell size or PIC and POC production rates measured in the same strains (Šupraha et al., 2015). This means that, when analyzed on the cellular and genotype level, coccolith volume cannot be used to indirectly infer cell size, carbon production rates or growth rates in this species. This observation confirms that intracellular coccolith production in Helicosphaera is constrained by genetics of biomineralization and defined by the dimensions of the organic template and coccolith vesicle, not by growth and other physiological rates (Henriksen et al., 2004; Young et al., 2004).

Similar to H. carteri, the mean CVI has no correlation to physiological rates or cell size in single genotypes of the heavily calcified species Coccolithus pelagicus and Coccolithus braarudii (Gerecht et al., 2014). Nevertheless, a strong positive relationship between coccolith dimensions and cell size of Coccolithus spp. is well established based on fossilized coccospheres, representing multiple genotypes and thousands of generations on geological timescales (Gibbs et al., 2018; Henderiks, 2008). Likewise, a positive correlation of coccolith size and cell size in Coccolithus pelagicus is evident when different genotypes from field communities and culture experiments are pooled together (Gibbs et al., 2013).

Apparently conflicting results have been reported for the members of the Noelaerhabdaceae family. For example, Fritz (1999) found that coccolith size does not correlate with physiological rates or cell size within a single genotype of Emiliania huxleyi, while studies by Bolton et al. (2016) and McClelland et al. (2016) found significant correlation between coccolith thickness, aspect ratio and the PIC : POC ratio within Noelaerhabdaceae, based on multiple strains from two closely related taxa (Emiliania and Gephyrocapsa). These observations support our conclusion that correlations between coccolith morphometry, cell size and physiology are weak or not identifiable within single strains (i.e., single genotypes) and that such relationships can only be established when considered across different genotypes or closely related species and lineages (Fig. 6).

Figure 6 Cross-plots of physiological and coccolith morphology parameters as measured in different strains of three modern coccolithophore species. (a) Cell volume (ln µm³) vs. coccolith volume (ln µm³). (b) Coccolith volume vs. aspect ratio (natural logarithmic scale). (c) Cell volume vs. PIC : POC ratio (natural logarithmic scale). (d) Coccolith aspect ratio vs. PIC : POC ratio. Error bars in (a), (b) and (d) represent ±1 SD of coccolith measurements in replicate cultures. Data sources: Šupraha et al. (2015; H. carteri) and McClelland et al. (2016; G. oceanica and E. huxleyi).

Across the modern representatives of the Helicosphaera genus, coccolith dimensions and cell size are correlated. For example, the smallest coccolith and coccosphere sizes are found in H. pavimentum (coccolith 4–6 µm, coccosphere 9–16 µm) and in H. hyalina (coccolith 5–8 µm, coccosphere 11–16 µm), whereas the largest coccolith and coccosphere sizes are found in H. carteri (coccolith 7–12 µm, coccosphere 15–25 µm) and H. wallichii (coccolith 7–12 µm, coccosphere 19–27 µm; http://www.mikrotax.org/Nannotax3/index.php?dir=Coccolithophores, last access: 4 May 2020; Geisen et al., 2004). It can therefore be argued with confidence that these traits followed the same principle in extinct fossil representatives, meaning that H. sellii had smaller coccolith and cell size than H. carteri, while H. granulata had the largest cell size of the three species. Thus, the variations in coccolith dimensions observed in the fossil record can be interpreted as an indicator of changing cell size both at the species and at the genus level.

4.2 Phenotypic evolution during the past 15 Myr – selection for smaller cells

The long-term phenotypic trends in fossil Helicosphaera indicate a selection for smaller cells in the South Atlantic Ocean and the equatorial Indian Ocean over the past 15 Myr. Our dataset supports that this selection for smaller cells occurred on various scales: (i) the (morpho)species level and (ii) the community level, as well as (iii) the evolutionary level through the extinction of larger species and speciation of smaller species.

At the (morpho)species level, this includes selection for smaller-celled phenotypes within each of the three most dominant species (H. carteri, H. sellii and H. granulata; Fig. 5). Along with the species-level trends, the selection for smaller cell size is also expressed in the relative contribution of different-sized species and their speciation patterns. This is evident in the gradual decrease in relative abundance and Pliocene extinction of the largest species H. granulata, whose dominance was gradually replaced by the rise in H. carteri and H. sellii. Subsequently, the extinction of H. sellii in the middle Pleistocene was followed by the rise in an even smaller extant species, H. hyalina. The adaptations at the species level and ecological trends on the community level, as well as the evolutionary processes (speciation and extinction) on the genus level, all followed the same pattern in favoring smaller-celled genotypes and species within Helicosphaera, a trend that is consistent between the two sites.

This selection for smaller morphotypes appears to be a common long-term adaptive response across multiple coccolithophore lineages. For example, the Calcidiscus leptoporus species complex (Calcidiscaceae) shows a decrease in the number of radial elements and coccolith diameter since the middle Miocene (Knappertsbusch, 2000; Fig. 7a), and selection for smaller species within the Noelaerhabdaceae family is well documented over the past 15 Myr across multiple deep-sea sites (Bolton et al., 2016; Hannisdal et al., 2012; Imai et al., 2015; Suchéras-Marx and Henderiks, 2014; Young, 1990; Fig. 7b). The consistent global response across coccolithophore lineages indicates that the phenotypic evolution of the group is driven by a common macroevolutionary driver (or drivers) that selected for smaller coccolith and thus smaller cell size since the middle Miocene.

Figure 7 Comparison of macroevolutionary phenotypic trends in three prominent coccolithophore lineages at temperate (blue) and tropical (red) sites, over the past 15 Myr. Mean coccolith size (length; µm) of (a) Helicosphaera (dots; this study) and Calcidiscus (open squares; Knappertsbusch, 2000) and (b) Reticulofenestra (triangles; this study in blue; Bolton et al., 2016, in red). Error bars are ±1 SD. (c) Mean aspect ratio (AR) in fossil specimens of Helicosphaera (blue and red dots; as in Fig. 4b) compared to AR estimates in Reticulofenestra spp. (grey triangles; grouped within size ranges of 2–3, 3–4, and 4–5 µm; Bolton et al., 2016), Emiliania and Gephyrocapsa (black symbols; McClelland et al., 2016). The light-blue shading encompasses the 1 SD variation in Helicosphaera over the entire time interval, whereas the long-term average AR would shift upwards if corrected for the systematic underestimation of thickness, as discussed in Sect. 2.3.

4.3 What drives the phenotypic evolution of Helicosphaera spp.?

The phenotypic expression observed at any point in the fossil record is a mean phenotype of thousands of generations of coccolithophore cells sedimented from the photic zone. It is a result of selection by global and local environmental drivers (e.g., hydrography, temperature, insolation, nutrients, CO₂ levels) as well as evolutionary (extinction, speciation) and ecological processes combined with inter- and intraspecific interactions such as competition and predation (Falkowski and Oliver, 2007; Finkel et al., 2007). Keeping in mind this environmental and biological complexity behind the final product – a coccolith with its phenotypic features – what does the selection for smaller cell size tell us about the possible global and regional drivers of phenotypic evolution in Helicosphaera and in coccolithophores in general?

The past 15 million years of Earth's geological history has included major environmental and paleoceanographic changes. Most importantly, this period saw the transition from the warm and high-CO₂ world of the Miocene climatic optimum (17–15 Ma) to the cooler, low-CO₂ world of the Pleistocene with its alternating ice ages and interglacial periods (e.g. Super et al., 2018; Zachos et al., 2001; Zhang et al., 2013). The global decrease in pCO₂ since the middle Miocene and associated carbon limitation of photosynthesis and calcification has been suggested as the main global driver of coccolithophore phenotypic evolution (Bolton et al., 2016). Carbon limitation during the late Miocene may have triggered a physiological adaptive response in coccolithophores around 7–5 Myr ago in the form of an active allocation of carbon from calcification to photosynthesis, as inferred from coccolith vital effects (Bolton and Stoll, 2013). A number of studies, spanning beyond the interval addressed here, have shown that carbon limitation especially affects larger, heavily calcified taxa, due to their low surface-to-volume (SA : V) ratio and high calcification demand for carbon (Hannisdal et al., 2012; Henderiks and Pagani, 2007; McClelland et al., 2016; Pagani et al., 2011). Being among the largest and most calcified species, members of the genus Helicosphaera would therefore be strongly affected by decreasing CO₂ levels during the middle to late Miocene, leading to progressive selection for smaller morphotypes with a higher SA : V ratio and lower carbon demand for calcification. However, a reversal in phenotypic trends is observed during the late Pliocene and early Pleistocene, with an increase in size despite a continuous cooling during the same period. Bolton et al. (2016) made the same observation in the Noelaerhabdaceae lineage and argued that such reversal in degree of calcification could be related to changes in ocean alkalinity which alleviated carbon limitation. In any case, the reversal in phenotypic evolution clearly indicates that other environmental drivers and adaptive responses were at play in the most recent time intervals.

While the global decrease in coccolith and cell size across coccolithophore lineages is likely related to decreasing pCO₂, the site-to-site contrasts observed in this study could be related to other environmental drivers which commonly shape coccolithophore communities such as nutrient availability and temperature. Experiments on H. carteri have shown that this species can strongly reduce its nutrient requirements, likely as an adaptation to growth in oligotrophic environments (Šupraha et al., 2015). This comes at the price of slower growth rates, on average larger cell size and slower calcification rates but conserved rates of photosynthesis, which is prioritized in respect to calcification under nutrient limitation. Similar short-term response patterns, where limited resources are allocated to maintain photosynthesis at the expense of calcification, were observed in semicontinuous cultures of C. pelagicus (Gerecht et al., 2015). On the other hand, E. huxleyi is much less sensitive to phosphorus limitation (Gerecht et al., 2018; Oviedo et al., 2014) due to its small size (and SA : V ratio), low nutrient requirements and high activity of its alkaline phosphatase (Riegman et al., 2000). These experimental observations show that nutrient availability can have a big impact on coccolithophores and that it likely exerts strong selective pressure in the natural environment, forcing the cells to channel limiting nutrients from calcification to photosynthesis to maintain their division rates or to develop physiological adaptations. Even though physiological low-nutrient adaptation was not reflected in the coccolith phenotype of single strains of H. carteri and C. pelagicus, it is plausible to expect that prolonged nutrient limitation would have strong effects on size selection over longer, evolutionary timescales (Falkowski and Oliver, 2007). In this sense, lower nutrient availability in the oligotrophic Indian Ocean as opposed to in the more productive South Atlantic (inferred from modern data as well as our % Fp-based PP estimates and nannofossil community composition) can account for the observed general offset in size between the two sites. In other words, lower nutrient concentrations coupled with carbon limitation could have channeled the Indian Ocean genotypes towards even smaller cell size compared to their South Atlantic counterparts. This process was likely supported and even exacerbated by consistently higher temperatures in the tropics (27±1 ^∘C, long-term average ±1 SD; Herbert et al., 2016) compared to in the South Atlantic and at higher latitudes (Fig. 2b). Cooler sea surface temperatures in the South Atlantic would have favored the solubility of CO₂ in the ocean, and coastal upwelling off SW Africa would have supplied nutrients to the photic zone, thus alleviating the effects of resource limitation and allowing for on average larger cell size at the South Atlantic site compared to at the Indian Ocean site.

An important adaptive trait to consider in the context of resource limitation is the heteromorphic, haplodiploid life cycle that is prevalent among coccolithophores including Helicosphaera (Young et al., 2005). The haploid life phase is generally significantly smaller and less calcified than the diploid life phase and is considered an adaptation to seasonal resource limitation (D'Amario et al., 2017; Šupraha et al., 2016). In addition, both life-cycle phases of Helicosphaera are flagellated and have a haptonema, which is a strong indication that they are putative mixotrophs (as demonstrated for Coccolithus pelagicus; Houdan et al., 2006) and can relieve nutrient limitation by ingesting bacteria. However, due to their sensitivity to dissolution, holococcoliths of Helicosphaera are scarce or near absent from the fossil record. It is thus difficult to assess how important the life-cycle dynamics of this genus were during its evolutionary history, but given the prevalence of this adaptive strategy in the modern ocean, it is safe to assume that it was also an important mechanism for overcoming resource limitation in the past.

4.4 Biogeochemical implications of phenotypic evolution

Despite the evolutionary selection for smaller coccolith and cell size in Helicosphaera since the middle Miocene, the trend in AR was less pronounced, indicating relatively consistent biogeochemical output (i.e., similar PIC : POC ratio; Fig. 8) of the genus. The observed AR trend suggests a slight increase in PIC : POC over time, which could mean that on the cellular level, the adaptive strategy of decreasing coccolith and cell size compensated for resource limitation and resulted in conserved physiological performance. Moreover, the regional trends in size selection were not reflected in regional patterns of AR, which were largely overlapping between both sites throughout the studied interval. Evidently, any global or regional patterns in phenotypic evolution, speciation or adaptation (either in size or in physiology) ultimately led to very similar biogeochemical output among the populations in different regions. A significant divergence in AR between the sites happened in the late Pleistocene, with an increase in AR (and thus PIC : POC) at Site 707. This coincided with an increase in PP estimate as well as coccolith size (Fig. 4), suggesting that locally weaker resource limitations could have led to increasing calcification rates and some degree of regional trends was possible despite overall globally conserved biogeochemical output.

Figure 8 Helicosphaera carteri (HC) coccolith dimensions in culture (black squares; mean ±1 SD) and fossil samples (grey circles; individual measurements at DSDP Site 525, N=1156) overlaid onto empirical relationship between coccolith size, coccolith thickness and PIC : POC (contours) based on culture data of E. huxleyi (EH) and G. oceanica (GO; McClelland et al., 2016).

However, when considered in a wider phylogenetic context and compared with AR patterns observed in the Noelaerhabdaceae (Fig. 7c), the AR trends observed in Helicosphaera, though statistically significant, appear biologically meaningless. During the same time interval, the AR and thus the biogeochemical output of Noelaerhabdaceae exhibits remarkable variability due to a wide range of possible PIC : POC ratios (i.e., degrees of calcification), with an overall decrease between the middle Miocene and Pleistocene. This phenotypic variation is likely related to the overall higher taxonomic and morphological diversity of the Noelaerhabdaceae, and it points to a high adaptive potential of this group that ensured its evolutionary success and global dominance during the Neogene, not just within coccolithophore communities (Figs. S3 and S4) but in marine phytoplankton in general. Our data show that, despite the global trends of size decreases across coccolithophore lineages, phenotypic evolution in different groups can lead to highly different biogeochemical outputs (Fig. 7c). A prime example of this is the comparison between the Helicosphaera genus, with its highly conserved morphology, little phenotypic innovation, low abundances and relatively stable biogeochemical output, and the Noelaerhabdaceae, with the potential for interspecies hybridization (Bendif et al., 2015) and rapid diversification (Bendif et al., 2019), global dominance in marine phytoplankton communities (Suchéras-Marx and Henderiks, 2014; Fig. S3), and a range of possible biogeochemical outputs (Figs. 6 and 8). Clearly, there are various adaptive strategies to the same global drivers among coccolithophore lineages and a wide range of biogeochemical consequences of their phenotypic evolution.

One possible explanation for different adaptive strategies can be related to the distinction between obligate calcifiers, including Coccolithus and Helicosphaera, and nonobligate calcifiers, such as Emiliania huxleyi (Durak et al., 2016; Walker et al., 2018). The phenotypic plasticity in obligate calcifiers could be more restricted when it comes to morphological innovation and biogeochemical output compared to nonobligate calcifiers, who can regulate their calcification rates (and thus PIC : POC ratio) to zero and overcome carbon limitation or other environmental stressors such as ocean acidification (Gafar et al., 2019). Obligate calcifiers, such as Coccolithus, Calcidiscus and Helicosphaera, combined were relatively more abundant than reticulofenestrids until ∼7 Ma when their fluxes decreased significantly at both sites (Fig. S3). This supports the idea that obligate calcifiers have been evolutionarily and ecologically less successful compared to nonobligate calcifier lineages during the global shift to the icehouse world of the late Miocene to Pleistocene, with each of the groups adapting within the constraints of their phenotype.

5 Conclusions

Inferring physiological adaptive strategies from coccolith morphometry in Helicosphaera and other coccolithophore taxa proves to be challenging if not impossible when considered within single genotypes. Valid correlations between morphology and physiology emerge only when more genotypes or closely related taxa are pooled together. Accordingly, short-term physiological responses in, e.g., cell size or PIC : POC ratios observed in single-genotype culture experiments do not necessarily translate into long-term phenotypic evolution observed in the fossil record, which is a result of far more complex and long-term selection by multiple biotic and abiotic drivers.

On multimillion-year timescales, environmental changes (directly or indirectly related to long-term cooling and a decrease in atmospheric CO₂ levels) caused selective pressures which resulted in a common decrease in cell size across different coccolithophore lineages, as herein shown for Helicosphaera (Zygodiscales order), Calcidiscus (Coccolithales order), and the ancestral lineage of Emiliania and Gephyrocapsa (Noelaerhabdaceae family; Isochrysidales order). The biogeochemical impact of the coccolithophores is defined by lineage-specific physiology and adaptive strategies, regionally distinct community compositions, and local ecophysiological adaptations. Over the past 15 million years, our data show that Helicosphaera displayed less phenotypic plasticity and a more conservative biogeochemical output compared to the Noelaerhabdaceae, which had (and still have) much higher variability in their degree of calcification. This supports the distinction between obligate and nonobligate calcifiers within coccolithophores, also on a macroevolutionary scale.