Haplo-diplontic life cycle expands coccolithophore niche

Coccolithophores are globally important marine calcifying phytoplankton that utilize a haplo-diplontic life cycle. The haplo-diplontic life cycle allows coccolithophores to divide in both life cycle phases and potentially expands coccolithophore niche volume. Research has, however, to date largely overlooked the life cycle of coccolithophores and has instead focused on the diploid life cycle phase of coccolithophores. Through the synthesis and analysis of global scanning electron microscopy (SEM) coccolithophore abundance data (n=2534), we find that calcified haploid coccolithophores generally constitute a minor component of the total coccolithophore abundance (≈ 2 %–15 % depending on season). However, using case studies in the Atlantic Ocean and Mediterranean Sea, we show that, depending on environmental conditions, calcifying haploid coccolithophores can be significant contributors to the coccolithophore standing stock (up to ≈30 %). Furthermore, using hypervolumes to quantify the niche of coccolithophores, we illustrate that the haploid and diploid life cycle phases inhabit contrasting niches and that on average this allows coccolithophores to expand their niche by ≈18.8 %, with a range of 3 %–76 % for individual species.

Our results highlight that future coccolithophore research should consider both life cycle stages, as omission of the haploid life cycle phase in current research limits our understanding of coccolithophore ecology. Our results furthermore suggest a different response to nutrient limitation and stratification, which may be of relevance for further climate scenarios.

Our compilation highlights the spatial and temporal sparsity of SEM measurements and the need for new molecular techniques to identify uncalcified haploid coccolithophores. Our work also emphasizes the need for further work on the carbonate chemistry niche of the coccolithophore life cycle.

1 Introduction

Coccolithophores are marine phytoplankton that produce calcium carbonate platelets, called “coccoliths”, which can be seen from space when coccolithophores bloom. Coccoliths eventually rain down into the ocean interior or serve as ballast as they are incorporated into faecal pellets and aggregates, which drives the carbonate pump and enhances the organic carbon pump by increasing organic carbon export rates to the deep sea (Klaas and Archer, 2002; Zeebe, 2012). Through the production of coccoliths, coccolithophores produce ≈1.5 Pg of inorganic carbon per year (Hopkins and Balch, 2018; Krumhardt et al., 2019) and subsequently account for 30 % to 90 % of carbonate in sediments (Broecker and Clark, 2009), highlighting the importance of coccolithophores in calcium carbonate burial.

In addition to the carbonate pump, coccolithophores contribute to the organic carbon pump, accounting for 1 %–40 % of marine primary production depending on habitat (Poulton et al., 2007, 2013). Because of involvement in the ocean carbon pumps and food web, coccolithophores thus play an important role in the ocean on regional to global spatial scales and seasonal to geological timescales.

Much focus has been put on understanding coccolithophore ecology and physiology, such as the function of calcification (Young, 1994; Monteiro et al., 2016; Xu et al., 2016), their diversity (Aubry, 2009; Young et al., 2003), and the factors controlling their calcification (Zondervan, 2007; Taylor et al., 2017) and competitiveness (Margalef, 1978; Krumhardt et al., 2017). However, one factor that significantly impacts coccolithophore calcite production and potentially their global success has been given little attention: their distinctive life cycle.

The life cycle of an organism is defined by the number of chromosome sets (the “ploidy level”) of the cell when asexual reproduction (“mitosis”) occurs. If mitosis occurs when the cell has one set of chromosomes (a haploid cell) the life cycle is called “haplontic” (Fig. 1a), while if mitosis occurs when the cell has two sets of chromosomes (a diploid cell) the life cycle is called “diplontic” (Fig. 1b). A few organisms can divide in both the haploid and diploid phase. Such a life cycle is called “haplo-diplontic” (Fig. 1c). Coccolithophores utilize the latter life cycle strategy – which is in contrast to dinoflagellates and diatoms, which tend to be either haplontic or diplontic and as such can only divide in either the haploid or diploid life cycle phase (Von Dassow and Montresor, 2011).

Figure 1 Life cycle strategies of phytoplankton: (a) dinoflagellates tend to utilize a haplontic life cycle, (b) diatoms tend to utilize a diplontic life cycle, and (c) coccolithophores tend to utilize a haplo-diplontic life cycle. Note that not all coccolithophores calcify in their haploid phase.

The haploid and diploid life cycle phases of coccolithophores can vary significantly in terms of coccolith structure, size, and morphology; cell size; and degree of calcification (Fig. 2). The diploid life cycle phases tend to be more heavily calcified than the haploid life cycle phases, which tend to be more lightly or non-calcified (Cros et al., 2000; Daniels et al., 2016; Fiorini et al., 2011a, b). This difference in cell calcium carbonate content (particulate inorganic carbon, PIC), cell organic carbon content (particulate organic carbon, POC) and the ratio thereof (the PIC : POC ratio) between the two life cycle phases means that the two phases potentially have contrasting impacts on the carbonate pump.

Figure 2 Coccosphere diversity of common coccolithophores (haploid cells are coloured in blue and diploid cells in red): (a) polycrater haploid morphology, (b) ceratolith haploid morphology, and (c–i) holococcolith haploid morphology. Note that in some instances multiple haploid phases are associated with one diploid phase (e.g. Syracosphaera mediterranea and Helicosphaera carteri), which may be due to cryptic speciation (Geisen et al., 2002). Furthermore, some species (e.g. E. huxleyi) do not calcify in their haploid phase and are thus not pictured. Images reproduced with permission from Young et al. (2020) (b–d, i) and Šupraha et al. (2016) (a, e–h). Panels (b)–(d) and (i) (HOL) were created by Jeremy Young, (i) (HET) was created by Marie-Helene Kawachi, and (a) and (e)–(h) were created by Luka Šupraha.

Although coccolithophore morphology is highly diverse, the diploid phases of coccolithophores primarily utilize heterococcolithophore morphology (with some exceptions, i.e. Braarudosphaera bigelowii), while the haploid life cycle phases can broadly be classified into four morphologies: polycrater (Fig. 2a), ceratolith (Fig. 2b), holococcolith (Fig. 2c–i), and unmineralized (not pictured) (Frada et al., 2018). Of these four haploid morphologies, the holococcolithophore morphology – which is defined by rhomboid calcite structures that constitute the coccoliths – is the most frequently utilized (Frada et al., 2018). Eight coccolithophore clades utilize holococcoliths, while four clades utilize an unmineralized haploid morphology, one clade utilizes a ceratolith morphology, and one clade utilizes ceratolith morphology, while for five clades the haploid morphology is currently unknown (Frada et al., 2018).

Coccolith and coccosphere morphology, cell and coccosphere size, and the degree of calcification influence coccolithophore ecology (Young, 1994). We can thus expect that the haploid and diploid life cycle phases of coccolithophores can have contrasting ecological preferences, which might allow a coccolithophore species to occupy multiple niches (Houdan et al., 2006; Frada et al., 2012; Cros and Estrada, 2013; Godrijan et al., 2018; Frada et al., 2018). This ability to occupy multiple niches should expand the total niche coccolithophore can inhabit, a potential advantage for haplo-diplontic organisms in variable environments (Mable and Otto, 1998). This is an idea that is supported by genetic models (Hughes and Otto, 1999; Rescan et al., 2015).

While niche differentiation has been widely observed for haplo-diplontic seaweeds (Couceiro et al., 2015; Guillemin et al., 2013; Lees et al., 2018; Lubchenco and Cubit, 1980) and coccolithophores (Houdan et al., 2006; Cros and Estrada, 2013; Godrijan et al., 2018; Frada et al., 2018), to date no research has quantitatively investigated the extent of niche overlap and niche expansion for haplo-diplontic algae. For coccolithophores this is because research has primarily focused on the diploid life phases, and relatively little is known in regard to the haploid life phase (Taylor et al., 2017; Frada et al., 2018). This is in part due to a research focus on the globally ubiquitous Emiliania huxleyi, which utilizes an unmineralized haploid morphology that cannot be readily identified with conventional light or scanning electron microscopy (Frada et al., 2008).

With the aim of understanding how haploid coccolithophores contribute to coccolithophore success, we quantify the niche overlap and niche expansion between haploid and diploid life stages of coccolithophores for the first time.

To do so, we compile global coccolithophore abundance observations of coccolithophores using all available scanning electron microscopy (SEM) measurements and where appropriate corresponding environmental measurements (temperature, salinity, dissolved inorganic nitrogen, phosphate, and silicate). Although our focus is on holococcolith-forming clades rather than E. huxleyi, holococcolith-forming clades include ecologically relevant species such as Helicosphaera carteri (Fig. 2e), Coccolithus pelagicus, and Calcidiscus leptoporus (Fig. 2h), which contribute more to the CaCO₃ flux to the deep ocean than E. huxleyi due to their larger coccolith and coccosphere size (Ziveri et al., 2007; Rigual Hernández et al., 2019).

In addition to niche overlap and niche expansion, we investigate the data set to identify ecological preferences of holococcolith-forming species, providing an updated picture on their global distribution, relative abundance, niche, and environmental controls. This work provides key information to better understand how the haplo-diplontic life cycle contributes to coccolithophore success.

2 Methods

2.1 Metadata compilation

Coccolithophore abundance measurements were compiled from 36 studies, constituting 2534 measurements and representing all major oceans (Table 1). These studies utilized scanning electron microscopy (SEM) to enumerate or further identify coccolithophores rather than solely relying on the more commonly utilized light or cross-polarized microscopy which under-represents coccolithophore biodiversity (Godrijan et al., 2018), especially in the case of holococcolithophores (Bollmann et al., 2002; Cerino et al., 2017). We used this data set to investigate global and vertical distribution patterns of haploid and diploid coccolithophore life cycle phases, specifically focusing on holococcolith forming species. Since abundance data were manually compiled, our data set is not exhaustive. For instance, some SEM studies, such as those by Okada and Honjo (1973), Honjo and Okada (1974), and Reid (1980), are not included in this data set since the data were not retrievable from the original publications.

Andruleit et al. (2003)Andruleit (2005)Andruleit (2007)Boeckel and Baumann (2008)Baumann et al. (2008)Cerino et al. (2017)Charalampopoulou et al. (2011)Charalampopoulou et al. (2016)Cepek (1996)Cros and Estrada (2013)D'Amario et al. (2017)Daniels et al. (2016)Dimiza et al. (2008)Dimiza et al. (2015)Eynaud et al. (1999)Giraudeau et al. (2016)Godrijan et al. (2018)Guerreiro et al. (2013)Guptha et al. (1995)Haidar and Thierstein (2001)Karatsolis et al. (2017)Kinkel et al. (2000)Luan et al. (2016)Malinverno (2003)Malinverno et al. (2015)Patil et al. (2017)Poulton et al. (2017)Saavedra-Pellitero et al. (2014)Schiebel et al. (2011)Schiebel et al. (2004)Smith et al. (2017)Šupraha et al. (2016)Takahashi and Okada (2000)Triantaphyllou et al. (2018)Silver (2009)

Table 1 Overview of metadata. The following abbrevations are used within this table: pLM stands for polarized light microscopy, LM stands for light microscopy, and SEM stands for scanning electron microscopy.

In addition to the global data set, we further investigated three case studies in order to better understand specific drivers and differences between the life cycle phases: the Atlantic Meridional Transect (AMT), representative of mid-oligotrophic open-ocean ecosystems; the long-term time series at Bermuda (BATS); and two time series in a mesotrophic coastal ecosystem in the Adriatic Sea (the “Mediterranean data set”). For the AMT study, we considered observations from four cruises, specifically AMT-12 (May–June 2003), AMT-14 (April–June 2004), AMT-15 (September–October 2004), and AMT-17 (October–November 2005), which had been previously published by Poulton et al. (2017). For the BATS station we considered data published by Haidar and Thierstein (2001), which consists of approximately monthly observations between January 1991 and January 1994. For the Mediterranean study, we combine two time series in the Adriatic Sea by Godrijan et al. (2018) and Cerino et al. (2017) taken between September 2008 to December 2009 and May 2011 to February 2013 at the RV-001 and C1-LTER stations, respectively.

For the BATS, AMT, and Mediterranean case studies, we additionally compiled temperature; salinity; and concentrations of dissolved inorganic nitrogen (DIN; nitrite + nitrate), phosphate, and silicate. For the AMT studies, environmental variables were acquired from the British Oceanographic Data Centre (BODC). For the BATS studies, environmental variables were acquired from the Bermuda Institute of Ocean Sciences (BIOS). For the Mediterranean study, day length was calculated using the MIT Skyfield package in Python. Other environmental variables such as turbulence, irradiance, and pH might also impact coccolithophore distribution patterns, but we have not included them in our compilation because they are not available for all presented case studies.

All data were acquired from supplementary data, online databases, or by contacting the authors directly if neither of these methods were available. The data were manually checked for synonyms or misspellings of species names, and where appropriate cell abundances were converted to cells L⁻¹. All species (or genera if not identified to a species level) were labelled as either heterococcolithophore, holococcolithophore, or “other”, which includes polycrater, nanoliths, and unidentified species. For these categorizations we followed definitions from Cros and Fortuño (2002).

The species and environmental data were compiled in Python and subsequently analysed in R (R Core Team, 2019). For all analyses we only considered samples within the top 200 m of the water column. To reduce the effects of seasonality, we binned the data into four main seasons, defined as December–February, March–May, June–August, September–November. We also calculated the mean of the observed abundances on a global scale and regional scale and estimated the highest observed abundances (the “maximum abundance”) for both heterococcolithophores and holococcolithophores and for each season. For the mean abundance calculations the mean was calculated for each sample and then averaged. Finally, we tested the count data for a normal distribution using a Shapiro–Wilk test for each region and the global data set. Where the count distribution was found to be normal (all data), a 95 % confidence interval was calculated.

Sampling bias and cover of data set

Our compilation contains sampling bias and is spatially and temporally incomplete. Temporally, there is bias towards the months June–August and December–February (29.28 % and 30.59 % of samples, respectively), with fewer samples in the inter-seasons. This temporal bias results from generally higher sampling effort in the Arctic Circle in June–August (8.43 % of samples) and the Southern Ocean in December–February (13.15 % of samples). Not coincidentally, this is when and where coccolithophore abundances are the highest (see results below). When excluding the Arctic Circle (June–August) and the Southern Ocean (December–February), the data set is temporally relatively evenly distributed (28.20 % March–May, 26.58 % June–August, 22.13 % September–November, 22.23 % December–February).

Spatially, there is higher sampling in the Atlantic Ocean, Mediterranean Sea, Arctic Circle, and Southern Ocean. In terms of spatial cover, coverage is limited in the Pacific Ocean and data is lacking in the Southern Ocean between June–August and the Arctic Circle between Dec–May. However, previous studies note the low coccolithophore abundance in the tropical and subtropical Pacific Ocean (Okada and Honjo, 1973; Honjo and Okada, 1974; Reid, 1980), and the absence or low abundance of holococcolithophores in this region (Okada and Honjo, 1973; Honjo and Okada, 1974; Reid, 1980). The lack of data in the Southern Ocean and the Arctic Circle for specific months is due to the difficulty of sampling these regions in the winter as well as low coccolithophore abundance due to light limitation.

The incomplete data cover of our data set combined with the spatial and temporal bias means that the analysis presented here mainly serves as a first-order estimate of the relative heterococcolithophore and holococcolithophore abundance and distribution patterns. For more accurate estimates, additional sampling needs to be conducted.

For more absolute estimates, additional sampling will have to be conducted. Inter-annual variability and strong links between coincident climate variability and primary productivity (Behrenfeld et al., 2006), as well as inter-annual and mesoscale variability on local scales, will influence phytoplankton distribution patterns (Volpe et al., 2012), which makes estimating global abundances challenging.

2.2 Definition of pairs and HOLP index

Not all heterococcolithophore-forming coccolithophore species form holococcospheres. Thus, to better illustrate the proportion of haploid and diploid coccolithophore cells, we reported the ratio between heterococcospheres and holococcospheres of species that form holococcoliths in their haploid phase, which is commonly implemented (Cros and Estrada, 2013; Šupraha et al., 2016).

This ratio is referred to as the “HOLP index” and is defined by Cros and Estrada (2013) as follows:

$$\begin{array}{}\text{(1)}& \mathrm{\text{HOLP index}}=\mathrm{100}\cdot {\displaystyle \frac{\text{paired holococcolithophore abundance}}{\text{paired coccolithophore abundance}}}.\end{array}$$

Species included in the HOLP index follow the definitions of paired species as defined in Frada et al. (2018) (Table 2), which are confined to currently understood associations and are likely to change as our understanding holococcolith species continues to improve. We calculated the HOLP index on a global and regional level for studies that identified holococcolithophores to a species level, the AMT data set, and the Mediterranean data set. To calculate the mean HOLP index, the ratios were calculated for each sample and then averaged.

Table 2 Taxonomic units included in HOLP index.

Taxonomic units included in HOLP index. Note that in some instances multiple heterococcolithophores are associated with single holococcolithophores (e.g. S. pulchra and S. protrudens are both associated with S. pulchra HOL).

2.3 Environmental drivers

We quantified the environmental drivers of heterococcolithophore and holococcolithophore abundance and the HOLP index for the AMT and Mediterranean data sets using Spearman correlations. We calculated Spearman correlations for heterococcolithophores and holococcolithophores and the HOLP index relative to temperature; salinity; depth; and concentrations of DIN (nitrite + nitrate), phosphate, and silicate for the AMT data set. The same ordinal associations were calculated for the Mediterranean data set, but we considered day length instead of depth because only the top 30 m of the water column was sampled and seasonality is an important driver in this region. To focus on marine systems of coccolithophores, we only considered samples with salinities above 30 ppt. Samples missing any environmental variables were removed. Subsequently, the AMT data set included a total of 45 samples, and the Mediterranean data set included 100 samples. Spearman correlation was performed in R using the “cor.test” function from the “stats” package (R Core Team, 2019). We also visualized environmental drivers by plotting the distributions of cell concentrations and environmental parameters within the water column or within the first two axes of a principal component analysis (PCA) and then interpolating values using the multilevel B-spline approximation (MBA) algorithm described by Lee et al. (1997). Prior to conducting the PCA, samples with a Cook's distance greater than 4 times the sample size were removed. For the visualizations, we used the same environmental parameters and samples as for the Spearman correlations, except for the AMT data set where we plotted chlorophyll instead of depth – which allowed for visualization of the deep chlorophyll maximum (DCM). For the AMT data set, we plotted the abundance and environmental parameters as a function of latitude and depth. While for the Mediterranean data set the variables were plotted as a function of the first two axes of a PCA, which included temperature; salinity; day length; and concentrations of phosphate, DIN, and silicate. Two different strategies were used to visualize the AMT and Mediterranean data sets, as the AMT data set is spatial and the Mediterranean data set is temporal. The MBA interpolation was performed with the “mba.surf” function from the “MBA” R package (Finley et al., 2017), and the PCA was performed with the “dudi” function of the “ade4” package (Dray and Dufour, 2007). Cook's distances were calculated using the “lm” and “cooks.distance” functions provided in the “stats” R package (R Core Team, 2019).

2.4 Seasonality

To investigate seasonality we compared monthly heterococcolithophore and holococcolithophore abundance data to temporal variations of temperature; salinity; day length; and concentrations of phosphate, DIN, and silicate of the BATS and Mediterranean data sets.

2.5 Niche overlap and niche expansion

Distribution patterns of phytoplankton are influenced by multiple environmental drivers. These environmental drivers form a n-dimensional hyperspace within which hypervolumes can be defined based on where the phytoplankton occur. This hypervolume is considered to be the species niche (Hutchinson, 1957) and allows niche comparisons between multiple phytoplankton – in this instance the two life cycle phases of coccolithophores.

Although processing hypervolumes is challenging due to their high dimensionality, methods described by Blonder et al. (2014) allow hypervolume quantification and comparison (for further discussion see Blonder, 2018 and Mammola, 2019). Using this strategy we determine the niche overlap of heterococcolithophores and holococcolithophores in hyperspace using the Sørensen–Dice and Jaccard similarity metrics.

We furthermore calculate the “niche expansion” of the haplo-diplontic life cycle strategy, which we define here as the non-overlapping region of either phase within hyperspace. In other words,

$$\begin{array}{}\text{(2)}& \mathrm{NE}\left(A\right)={\displaystyle \frac{\left|A\right|-|A\cap B|}{|A\cup B|}},\end{array}$$

where NE(A) is a niche expansion of A, A is hypervolume A, B is hypervolume B, ∩ is the intersection between two hypervolumes, and ∪ is a union between two hypervolumes.

The niche metrics utilized in this study are illustrated in Fig. 3. Although we visualize the niche of each species using contours, in reality the niche metrics are calculated based on random points sampled from the inferred hypervolumes (Blonder et al., 2014).

Figure 3 Hypervolume metrics utilized in this study: (a) union, (b) intersection, (c) unique fraction A, (d) Jaccard similarity metric, (e) niche expansion, and (f) Sørensen–Dice similarity.

We calculated the Jaccard and Sørensen–Dice similarity metrics and niche expansion for the AMT, BATS, and Mediterranean Sea data sets. For the AMT data set, DIN showed high Pearson correlation to silicate (ρ=0.95, p<0.001) and phosphate (ρ=0.90, p<0.001). We thus only considered temperature, salinity, and the concentration of DIN in this region. Although no such correlation was observed for the Mediterranean data set, and weaker but significant relationships were observed in the BATS stations (ρ=0.74, p<0.001 for silicate and ρ=0.84, p<0.001 for phosphate), to make the niche metrics comparable in all regions the silicate and phosphate concentrations of the Mediterranean and BATS data sets were also excluded.

It is likely, however, that silicate and phosphate, as well as other parameters (such as irradiance, turbulence and carbonate chemistry), influence the niche of coccolithophores and thus the metrics calculated. Besides the influence of environmental parameter choice, results of the niche analysis will depend on what is considered a paired species. Although we use up-to-date definitions from Frada et al. (2018), these definitions are likely to change in the future. Finally, cryptic speciation (Geisen et al., 2002) and subsequently the pairing of multiple haploid holococcolith (HOL) phases to single diploid heterococcolith (HET) phases and vice versa complicate results.

The environmental data were normalized using z scores prior to analysis. Niche overlap and niche expansion were calculated only for species for which both life cycle phases were observed. In addition to calculating the niche expansion for individual species, we calculated an average niche expansion by taking the mean NE values of all individual species for both the haploid and diploid coccolithophore life cycle phases.

We used the “hypervolume” R package (Blonder and Harris, 2018) to conduct our niche overlap and niche expansion analysis. Gaussian kernel density estimation (R function “hypervolume_gaussian”) was used to construct the hypervolume, the overlap metrics were calculated with the “hypervolume_overlap_statistics” R function, and the volume and intersection of hyper volumes were calculated using the “get_volume” R function.

3 Results

3.1 Biogeography of coccolithophores

Within our compilation, heterococcolithophores showed global distribution, while holococcolithophores were noticeably absent at the ALOHA station in Hawaii and (with some exceptions) >50^∘ S in the Southern Ocean (Fig. 4 and Table 3).

Figure 4 (a–b) Global coccolithophore distribution and (c–h) latitudinal coccolithophore distribution: (a) Heterococcolithophores, (b) Holococcolithophores, (c) Heterococcolithophores, (d) E. huxleyi, (e) paired heterococcolithophores, (f) holococcolithophores, (g) paired holococcolithophores, and (h) HOLP index. For the latitudinal plots, the light shading is log-transformed distribution.

Table 3 Global heterococcolithophore and holococcolithophore abundance.

Values in parentheses are 95 % confidence intervals.

The highest maximum abundances of heterococcolithophores are observed at high latitudes within the Arctic Circle (>66^∘ N) ($\approx \mathrm{4.37}\times {\mathrm{10}}^{\mathrm{6}}$ cells L⁻¹ for June–August) and the Southern Ocean (>40 and <65^∘ S) ($\approx \mathrm{1.64}\times {\mathrm{10}}^{\mathrm{6}}$ cells L⁻¹ for December–February). Generally, maximum abundances above 1×10⁵ cells L⁻¹ were observed, except between September–November in the Indian Ocean (≈ 3.33×10⁴ cells L⁻¹), September–November and December–February in the Atlantic Ocean ($\approx \mathrm{5.40}\times {\mathrm{10}}^{\mathrm{4}}$ and $\approx \mathrm{9.78}\times {\mathrm{10}}^{\mathrm{4}}$ cells L⁻¹ respectively), and March–May in the Pacific Ocean ($\approx \mathrm{4.96}\times {\mathrm{10}}^{\mathrm{4}}$ cells L⁻¹).

The regions and periods with the highest mean heterococcolithophore abundance differ from the regions and periods with the highest maximum heterococcolithophore abundance. For example, the highest mean abundance is observed in the Indian Ocean during March–May ($\approx \mathrm{1.13}\times {\mathrm{10}}^{\mathrm{5}}$ ($\pm \mathrm{2.97}\times {\mathrm{10}}^{\mathrm{4}}$) cells L⁻¹), which is higher than the highest mean abundance in the Southern Ocean observed during March–May ($\approx \mathrm{1.17}\times {\mathrm{10}}^{\mathrm{5}}$ ($\pm \mathrm{2.88}\times {\mathrm{10}}^{\mathrm{4}}$) cells L⁻¹) and in the Arctic Circle during June–August ($\approx \mathrm{5.83}\times {\mathrm{10}}^{\mathrm{4}}$ ($\pm \mathrm{2.97}\times {\mathrm{10}}^{\mathrm{4}}$) cells L⁻¹).

Although holococcolithophores show low abundances in the high latitudes of the Southern Hemisphere, highest maximum holococcolithophore abundances are observed in the Arctic circle (>66^∘ N) during June–August ($\approx \mathrm{2.23}\times {\mathrm{10}}^{\mathrm{5}}$ cells L⁻¹). High maximum abundances are additionally observed in the Mediterranean Sea (September–November) ($\approx \mathrm{1.27}\times {\mathrm{10}}^{\mathrm{5}}$ cells L⁻¹).

The lowest maximum holococcolithophore abundance is observed in the Pacific Ocean during June–August (4.45×10² cells L⁻¹) and in the Arctic Circle during September–November (1.12×10³ cells L⁻¹).

On average, the Mediterranean Sea has the highest mean holococcolithophore abundance (between $\approx \mathrm{2.21}\times {\mathrm{10}}^{\mathrm{3}}$ and 9.42×10³ cells L⁻¹), followed by the Indian Ocean ($\approx \mathrm{1.41}\times {\mathrm{10}}^{\mathrm{3}}$–4.80×10³ cells L⁻¹). The lowest mean abundances are observed in the Pacific Ocean (4.9×10¹ ($\pm \mathrm{9.70}\times {\mathrm{10}}^{\mathrm{1}}$), Arctic Circle (September–November; 2.55×10² ($\pm \mathrm{2.71}\times {\mathrm{10}}^{\mathrm{2}}$), and Southern Ocean (December–February; $\approx \mathrm{3.24}\times {\mathrm{10}}^{\mathrm{2}}$ ($\pm \mathrm{2.06}\times {\mathrm{10}}^{\mathrm{2}}$) cells L⁻¹).

Depending on the season, holococcolithophore contribution to total coccolithophore abundance varies globally between 1.67 % (±0.37 %) in December–February and 16.16 % (±1.68 %) in June–August, with the highest contribution observed in the Mediterranean Sea in June–August (31.38 % ± 2.93 %) (Table 3). On an regional scale outside of the Mediterranean Sea, holococcolithophores contribute less than 8 % to the total coccolithophore abundances. However, the contribution of holococcolithophores to paired species is higher than when all heterococcolithophore and holococcolithophores are considered (Table 4), with a HOLP index between 5.65 % (±1.71 %) and 27.41 (±2.67 %) globally depending on season. The lowest HOLP indices were observed in the Atlantic Ocean in September–November (0.59±0.81) and December–February (0.47±0.65 %), and in the Southern Ocean in December–February (0.61±0.58 %). The highest HOLP index was observed in the Mediterranean Sea in June–August (39.03±3.23 %).

Table 4 Global HOLP index.

Mean HOLP indices grouped by season and location. Values in parentheses are the 95 % confidence interval.

3.2 Vertical distribution

In the global data set, heterococcolithophore abundance is evenly distributed with depth, while holococcolithophore abundance is highest in the top 50 m of the water column (Fig. 5).

Figure 5 Global depth distribution of heterococcolithophores and holococcolithophores: (a–d) total paired and unpaired heterococcolithophore and holococcolithophore abundance and (e–f) individual species abundances. Heterococcolithophores are plotted in red, and holococcolithophores are plotted in blue. Only the most abundant coccolithophore species are plotted individually. Error bars are the standard error.

For holococcolithophores the vertical distribution pattern is mainly driven by paired holococcolithophore species, which constituted ≈62.2 % of the total coccolithophore abundance. Two currently unpaired holococcolithophores also contribute to the depth distribution trend with Helladosphaera cornifera (for which the association has to be further confirmed), constituting ≈8.1 % of total holococcolithophore abundance, and Corisphaera gracilis (for which no pair has been described), constituting ≈3.6 % of total holococcolithophore abundance. Subsequently paired holococcolithophore abundances broadly followed the same patterns observed when all holococcolithophores were considered.

In comparison to holococcolithophores, depth distribution of heterococcolithophores was driven by unpaired species, in particular E. huxleyi, which constituted ≈59.2 % of total heterococcolithophore abundance, but was also driven by the presence of unpaired deep-water species such as Ophiaster formosus, Florisphaera profunda, Calciopappus caudatus, and Oolithotus antillarum. However, although paired heterococcolithophores only contributed ≈5.7 % to total heterococcolithophore abundance, the depth distribution trends of paired and total heterococcolithophores species were similar.

3.3 Environmental drivers of niche partitioning

To further understand the distribution patterns observed on a global basis and within the water column we investigated the environmental drivers of heterococcolithophore and holococcolithophore abundance in the Atlantic Ocean (with the AMT data set) and the Mediterranean Sea. For the Atlantic Ocean data set, the environmental drivers were considered in the context of their distribution within the water column, whereas for the Mediterranean the environmental drivers were considered within PCA “niche space”. These observed patterns were then further corroborated through Spearman analysis.

3.3.1 Atlantic Ocean

In the Atlantic Ocean both heterococcolithophores and holococcolithophores have their highest abundances in the top 50 m of the water column (Fig. 6). However, a noticeable difference between heterococcolithophore and holococcolithophore distribution (Fig. 6a and d respectively) is the absence of holococcolithophores below the deep chlorophyll maximum (DCM) (Fig. 6l). The DCM tends to occur at 1 %–10 % irradiance levels and is closely linked to the nutricline and thermocline (Poulton et al., 2006). The difference in depth distribution between heterococcolithophores and holococcolithophores and the absence of holococcolithophores below the DCM may therefore be influenced by a combination of light limitation, high nutrient concentrations, cold water temperatures at depth, or other factors not addressed in this study.

Figure 6 Depth distribution along AMT: (a) heterococcolithophore abundance, (b) paired heterococcolithophore abundance, (c) E. huxleyi abundance, (d) holococcolithophore abundance, (e) paired holococcolithophore abundance, (f) HOLP index, (g) temperature (^∘C), (h) salinity (ppt), (i) DIN (µM), (j) silicate (µM), (k) silicate, and (l) chlorophyll. Species abundances are plotted on log scale.

This suggests that heterococcolithophores might be better adapted to exploit such conditions. Although differences in sinking rates, which are conceivably higher in the more heavily calcified heterococcolithophores, could also factor into the difference in depth distribution between the two life cycle phases.

Table 5 PCA loadings of the Mediterranean Sea data set.

The first two axis of the PCA captured 53.94 % of variance. Data are from Cerino et al. (2017) and Godrijan et al. (2018).

The distribution of heterococcolithophores (Fig. 6a) is primarily driven by E. huxleyi (Fig. 6c), which constitutes ≈30 % of total heterococcolithophore abundance in the data set. When only paired heterococcolithophore species were considered (Fig. 6b), a more even distribution in subtropical and tropical regions is observed. Holococcolithophores and paired holococcolithophores showed roughly similar distribution patterns (Fig. 6d and e).

Within the upper water column, heterococcolithophores showed the highest abundance at higher latitudes (>35^∘ N and >30^∘ S), which is associated with a shallow mixed layer, lower salinity, and lower temperature, as well as increasing silicate concentrations in the Southern Hemisphere. Holococcolithophores meanwhile showed the highest abundances at both high latitudes and in the Atlantic subtropical gyres. The HOLP index (Fig. 6f) was highest within the Atlantic subtropical gyres, with a higher proportion of holococcolithophores in the northern subtropical gyre, which is associated with a shallower DCM relative to the southern subtropical gyre. This shallowing of the DCM on the AMT is, however, likely a seasonal signal as described by Poulton et al. (2006) and Poulton et al. (2017).

Spearman correlations (Table 6) suggests holococcolithophores are significantly (p<0.05) negatively correlated to phosphate, DIN, silicate, and depth and significantly positively correlated to temperature and salinity. Paired holococcolithophores and the HOLP index showed the same correlation trends as holococcolithophores.

Table 6 Spearman correlations for the AMT data set.

${}^{***}$ p<0.001. ${}^{**}$ p<0.01. ^* p<0.05. Significant correlations are highlighted in bold. Data were acquired from Poulton et al. (2017). The following abbreviations are used within this table: HET stands for heterococcolithophores, HET.P stands for paired heterococcolithophores, HOL stands for holococcolithophores, HOL.P stands for paired holococcolithophores, and HOLP stands for the HOLP index.

On the contrary, heterococcolithophores are only significantly and negatively correlated with depth and phosphate. While for paired heterococcolithophores significant negative correlations were observed with depth and silicate.

Thus heterococcolithophore and holococcolithophore abundance in the Atlantic Ocean seems primarily driven by the depth of the DCM both in terms of vertical and latitudinal distribution. The highest abundances of both heterococcolithophores and holococcolithophores are observed above the DCM, and heterococcolithophores are present below the DCM while holococcolithophores are not. In terms of latitude, the highest abundances of heterococcolithophores correspond to the shallow DCM depth that occurs in higher-latitude regions, and the highest abundances of holococcolithophores occur in subtropical regions with deep DCM depths.

3.3.2 Mediterranean Sea

For the Mediterranean Sea long-term time series, niche separation of heterococcolithophores and holococcolithophores within the PCA niche space (Fig. 7) is primarily driven by principal component 1 (PC1), which is positively associated with temperature and day length and negatively associated with salinity, DIN, silicate, and phosphate (see Table 7). Heterococcolithophores are most abundant at low PC1 values (i.e. the left quadrants of Fig. 7a), which correspond to low temperatures and short day lengths and high salinity and concentrations of DIN, silicate, and phosphate (see Table 5). Holococcolithophores are most abundant at high PC1 values (i.e. the right quadrants of Fig. 7b), which correspond to high temperatures and long day lengths and low salinity and concentrations of DIN, silicate, and phosphate.

Figure 7 Principal component analysis (PCA) of the RV-001 and LTER1 stations in the Mediterranean Sea. Abundance and environmental values were projected on the PCA post hoc and then interpolated: (a) heterococcolithophore abundance, (b) holococcolithophore abundance, (c) salinity, (d) temperature, (e) depth, (f) phosphate, (g) DIN, and (h) silicate. Data were acquired from Cerino et al. (2017) and Godrijan et al. (2018).

Table 7 Spearman correlations for the Mediterranean data set.

${}^{***}$ p<0.001. ${}^{**}$ p<0.01. ^* p<0.05. Significant correlations are highlighted in bold. Data were acquired from Godrijan et al. (2018) and Cerino et al. (2017).The following abbreviations are used within this table: HET stands for heterococcolithophores, HET.P stands for paired heterococcolithophores, HOL stands for holococcolithophores, HOL.P stands for paired holococcolithophores, and HOLP stands for the HOLP index.

The pattern observed in the PCA niche space should be interpreted with some caution because only a portion of the variance is captured (53 %) and the use of interpolation introduces additional uncertainties. Besides, the structure of the PCA depends highly on the number and type of variables included (Figs. S2–S4 in the Supplement), particularly when time is considered. However, the patterns presented in the PCA are also apparent in the Spearman correlations (see Table 6), which suggests that the PCA plots are a qualitatively good representation of the data.

The Spearman correlations indicate that heterococcolithophores are significantly negatively correlated to temperature and day length and significantly positively correlated to phosphate, DIN, silicate, and salinity. For paired heterococcolithophore species the only significant correlation observed was a positive correlation with silicate.

Holococcolithophores showed the opposite pattern to heterococcolithophores and are significantly positively correlated to day length and temperature and significantly negatively correlated to salinity, DIN, silicate, and phosphate. Paired holococcolithophores and the HOLP index showed significant positive correlation to temperature and day length, but no significant correlations with the other environmental variables were observed.

3.3.3 General environmental trends

Our statistical analysis shows that in both the Mediterranean Sea and Atlantic Ocean holococcolithophores are generally found in low-nutrient and warm environments and high light availability. However, an opposite trend was observed between the Atlantic Ocean and Mediterranean Sea in terms of correlation to salinity, with holococcolithophores positively correlated to salinity in the Atlantic Ocean and negatively correlated to salinity in the Mediterranean Sea. This difference in correlation to salinity may be explained by the different drivers of salinity in both regions. In the Atlantic Ocean, low salinity occurs at high latitudes, while high salinity corresponds to mid-ocean gyres due to higher evaporation in tropical and subtropical regions. In contrast, at the coastal site in the Mediterranean Sea, low salinity is strictly related to direct freshwater input and associated nutrients. As such salinity may be simply correlated to other environmental drivers, rather than be a driver itself.

In the Mediterranean Sea and the Atlantic Ocean significant negative correlations were observed between holococcolithophores and silicate. Although this correlation could be in part due to strong correlation between DIN and silicate (ρ=0.95) observed in the Atlantic Ocean, the reason for this is less clear in the Mediterranean Sea as no such correlation is observed. A physiological reason for the negative correlation to silicate could be different silicate requirements among different coccolithophore species. Durak et al. (2016) for instance found evidence of silicate requirement for the heterococcolith life cycle phases of S. apsteinii, C. coccolithus, and C. leptoporus but not for E. huxleyi or G. oceanica. Follow-up experiments have furthermore found holococcolith life cycle phases of C. coccolithus and C. leptoporus do not require silicate (Langer et al., 2021).

Statistically significant correlations were the same when all holococcolithophores, paired holococcolithophores, or the HOLP index was considered at both locations; however, fewer significant correlations were observed for paired holococcolithophores and the HOLP index.

The trend for heterococcolithophores is less clear when comparing the two sites: an opposite trend to holococcolithophores, i.e. high nutrients and low temperatures, is observed in the Mediterranean Sea but not in the Atlantic Ocean where many of the correlations were not significant and heterococcolithophore were negatively correlated to phosphate. This negative correlation to phosphate is potentially due to deeper sampling in the Atlantic Ocean combined with high phosphate concentrations in deep and light-limited waters skewing correlations, which highlights the need to consider sampling and DCM depth when comparing environmental correlation between studies. It may furthermore be due to the presence of mixotrophic or heterotrophic coccolithophores at depth in the Atlantic Ocean, which are not found in the shallow coastal waters of the Mediterranean Sea.

3.4 Niche overlap and niche expansion

We conducted niche similarity and niche expansion calculations on the AMT, BATS, and Mediterranean data sets to quantify niches in these regions. For niche overlap we considered the Jaccard overlap and Sørensen–Dice overlap metrics, which range from 0 to 1, with 1 signifying complete overlap. For niche expansion we considered the relative amount each life cycle contributed to the total niche volume. In the AMT the niche overlap of paired species was high for both the Jaccard overlap and Sørensen–Dice overlap metrics (0.84 and 0.91, respectively; see Table 8). However, for individual species the overlap metrics were highly variable, ranging from 0.11 to 0.74 and from 0.20 to 0.81 for the Jaccard overlap and Sørensen–Dice overlap metrics, respectively. The niche expansion was higher for heterococcolithophores than holococcolithophores when all paired species were considered (see Table 8) but was again highly variable for individual species. The holococcolithophore phase of C. mediterranea, S. bannockii, H. wallichii, and C. leptoporus for instance all contributed more to the total niche volume than their heterococcolithophore life cycle phase.

Table 8 Niche expansion (NE) and niche overlap.

Niche overlap and niche expansion metrics utilized in this study. For definitions, see Sect. 2.5 and Fig. 3.

For BATS, the niche overlap values are generally smaller than for the AMT, with a Jaccard overlap and Sørensen–Dice values of 0.60 and 0.66, respectively. The niche expansion of heterococcolithophores at BATS is larger compared to the AMT (0.49 versus 0.11 for BATS and AMT, respectively). The NE of holococcolithophores is similar for both stations (0.02 versus 0.05 for BATS and the AMT, respectively). S. anthos and S. pulchra are the only species for which coccolithophore life cycle pairs are observed at the BATS station. For these species, the NE of heterococcolithophore is similar to when all species were considered but is higher for holococcolithophores. In the Mediterranean Sea, the niche overlap and niche expansion values are more similar to the BATS data set than to the AMT data set.

Niche expansion of heterococcolithophores was also higher than holococcolithophores when all paired species were considered, but like in the Atlantic Ocean species-specific exceptions were observed. The holococcolithophore phase of C. mediterranea, S. histrica, S. strigilis, and C. leptoporus all contributed more to the total niche volume than their heterococcolithophore life cycle phase in this region. In the Mediterranean Sea the niche of S. molischii is of particular note, as no overlap between the two life cycle phases was observed, and the two unique components were of similar size (0.51 and 0.49 for heterococcolithophores and holococcolithophores, respectively).

Although quantitative interpretation of niche is difficult since niche will vary depending on the number of environmental axes included (Blonder et al., 2014), these results highlight that holococcolithophores contribute significantly to the niche volume of coccolithophores, in some instances contributing more to total niche volume than the heterococcolithophore phase. In this context C. pelagicus is particularly relevant as this species contributes significantly to the global carbonate flux (Ziveri et al., 2007; Rigual Hernández et al., 2019 and is one of the key calcifiers in the Arctic Ocean (Daniels et al., 2016).

These results additionally suggest that the niche expansion patterns of the coccolithophore life cycle are more similar between the BATS and Mediterranean Sea than BATS and the AMT. This suggest that seasonal variations play an important role in structuring the niche of coccolithophores, otherwise BATS and the AMT should be more alike due to more similar hydrographic conditions. This result highlights the value of time series for studying the ecology of the coccolithophore life cycle and raises the need for caution when comparing niche volumes of cruise data and time series.

3.5 Seasonality of coccolithophores

Heterococcolithophore and holococcolithophore abundance highly varies with season at both the BATS station in the Atlantic Ocean and the long-term stations in the Mediterranean Sea (Fig. 8). Both locations experience a peak of heterococcolithophores in the winter, followed by a peak of holococcolithophores at the end of spring and in early summer. In the Atlantic Ocean, the heterococcolithophore are present in high abundance for a longer period of time, overlapping with the spring peak in holococcolithophore abundance (Fig. 8b).

Figure 8 Seasonality of heterococcolithophores and holococcolithophores at the BATS station in Bermuda (left column) and the RV-001 and LTER-1 stations in the Mediterranean Sea (right column). Note that heterococcolithophores are most abundant in the winter, followed by a high abundance of holococcolithophores in the late spring and early summer: (a–b) heterococcolithophore abundance, (b–c) holococcolithophore abundance, (e, g) temperature, (f, h) chlorophyll, (i, k) DIN (nitrite + nitrate), (j, l) silicate.

At both locations the peak of the holococcolithophore bloom occurs in the spring and summer when water temperatures rise and the day length is longest, while heterococcolithophore abundance is highest in the winter when temperature is lowest and day length shortest. The seasonality of peak heterococcolithophore and holococcolithophore abundance may furthermore correspond to seasonal changes in mixed layer depth (MLD), as both the Atlantic Ocean and Mediterranean Sea experience increased mixing in the winter and higher stratification in the summer.

Clear seasonal patterns are observed for DIN and silicate concentrations at both locations. High holococcolithophore abundances are observed when silicate and DIN concentrations are low and vice versa for high heterococcolithophore abundances. This observed seasonal patterns at the BATS and Mediterranean time series are thus consistent with our Spearman correlations (see Tables 6 and 7 and our discussion above) and PCA niche space (see Fig. 7 and discussion above).

It is important to note that on a species level individual species do not exclusively follow the seasonal heterococcolithophore or holococcolithophore trends described above, as illustrated in detail previously (Cerino et al., 2017; Godrijan et al., 2018). For instance, for Syracosphaera molischii and Syracosphaera pulchra the holococcolith phase rather than heterococcolith phase is the dominant life cycle phase in these time series. Furthermore, the holococcolithophore phases of S. molischii, Syracosphaera histrica, Algirosphaera robusta, and Acanthoica quattrospina are observed in the winter – a period when total holococcolithophore abundance is lowest. Finally, on a individual level succession does not immediately follow the previous life cycle phase, with several months of absence observed between peak abundance for some species (Cerino et al., 2017; Godrijan et al., 2018).

This highlights that grouped heterococcolithophore and holococcolithophore abundances represent a generalization that might not always represent patterns observed for individual species. These differences from generally observed patterns could be due to variations in life strategy – such as mixotrophy, motility, and grazing susceptibility – independent of life cycle phase, suggesting that functional traits different from the life cycle phase may determine the niche these species inhabit.

4 Discussion

Our meta-analysis shows that holococcolithophores are a minor contributor to coccolithophore abundance in the modern ocean, contributing between ≈ 2 %–15 % to the total coccolithophore abundance and between ≈ 5 %–30 % of the total paired coccolithophore abundance depending on season. However, our analysis also shows that haploid cells play an important role in coccolithophore ecology, accounting for ≈19 % of their niche volume, with lesser or greater contributions depending on the species (3 %–76 %). Our analysis furthermore shows that if conditions are favourable (specifically increased stratification and reduced nutrient supply) holococcolithophores can be significant contributors to the coccolithophore standing stock (up to ≈30 %).

Although holococcolithophore contribution to calcium carbonate production is likely small due to their lower cellular CaCO₃ content – which is an order of magnitude lower than heterococcolithophores (Daniels et al., 2016; Fiorini et al., 2011a, b) – their role in the carbonate cycle in present, past, and future oceans could have other biogeochemical effects. A shift towards a higher proportion of holococcolithophore cells would result in lower global calcium carbonate production, which could subsequently result in lower CO₂ outgassing on short timescales. Furthermore, the ballasting effect of coccolithophores would be reduced if a shift towards more lightly calcified haploid cells occurred (Hoffmann et al., 2015), which would potentially reduce efficiency of the carbon pump by reducing sinking rates. Although how other factors such as shifts in carbonate chemistry impact holococcolithophore abundance are not clear, increased stratification and decreased nutrient supply are projected under the RCP 8.5 climate change scenario (Fu et al., 2016), which would favour holococcolithophores. This shift from diploid to haploid coccolithophores could, on the one hand, reduce CO₂ outgassing but would, on the other hand, additionally reduce ballasting and subsequently impact the carbon pump by reducing sinking rates.

In terms of the ecological niche – which is the environmental range a species inhabits – observations of heterococcolithophores and holococcolithophores in our meta-analysis broadly conform to the Margalef niche space model. This model was proposed by Margalef (1978) and posits that the distribution of phytoplankton functional groups relate broadly to turbulence, light, and nutrients. Although we do not explicitly represent the former, turbulence is implicitly represented in our analysis based on mixed-layer depth. Within the Margalef niche framework, we find that heterococcolithophores and holococcolithophores occupy an intermediate functional group located between diatoms and dinoflagellates (see Fig. 9), as proposed by Houdan et al. (2006) and Frada et al. (2018).

Figure 9 A modified version of Margalef's niche model (Margalef, 1978), as proposed by Houdan et al. (2006) and Frada et al. (2018). Note that we have added day length, which was proposed by Balch (2004) as a third axis for the Margalef niche model.

Diploid heterococcolithophores thus favour high-nutrient, more turbulent waters, whereas haploid holococcolithophores favour low-nutrient, more stratified waters.

Although the Margalef niche model certainly presents a simplification that is prone to exceptions both for diatoms (see Kemp and Villareal, 2018) and coccolithophores (the generalist E. huxleyi and deep-water species such as F. profunda are clear examples), the model broadly holds in our meta-analysis.

This ecological–environmental distinction of heterococcolithophores and holococcolithophores is observed in coccolithophore species distribution in our analysis in terms of geographical succession, depth distribution and seasonal trends. In the Southern Ocean and Atlantic Ocean a geographical shift from holococcolithophores to heterococcolithophores is observed as latitude and turbulence and nutrients increases, while in the Atlantic Ocean and the global data set a vertical shift is observed, with holococcolithophores absent or at low abundance in deep nutrient-rich waters. Finally, in the Mediterranean Sea a seasonal shift is observed, as heterococcolithophores are most abundant in well-mixed nutrient-rich winter months and holococcolithophores are most abundant in nutrient-poor stratified summer months.

However, some exceptions occur. For instance, in the AMT data set, although heterococcolithophores are more evenly distributed with depth, the maximum abundance of heterococcolithophores is in surface waters, and subsequently heterococcolithophores are negatively correlated to nutrients. Nonetheless, the relation to turbulence holds: heterococcolithophore abundance is highest in well-mixed high-latitude waters, and holococcolithophore abundance is highest in stratified subtropical regions. Finally, many species-specific exceptions occur. We highlight examples on a seasonal scale in our Mediterranean data set discussion (see Sect. 3.6), but exceptions were also noted along the AMT (see discussion in Poulton et al., 2017) and in other Mediterranean studies (Šupraha et al., 2016; D'Amario et al., 2017; Skejić et al., 2018). This means that caution should be used when considering the niche model for individual species.

4.1 Niche overlap and expansion

Our study showed that the niche volume of coccolithophores is larger when holococcolithophores are included in coccolithophore niche volume. This tells us two things: first, studies focused solely on heterococcolithophores are underestimating coccolithophore habitat and thus inaccurately represent the coccolithophore functional group in modelling and physiological studies. This means that we might be underestimating their ability to compete with other phytoplankton and the range of environmental conditions that they can tolerate. Secondly, we underestimate coccolithophore primary productivity and calcite standing stock by not including accurate assessments of their abundance.

This might be of particular relevance for E. huxleyi, the diploid phase of which has been of particular research focus due to its high abundance (approx 59.2 % in our compilation). Although our meta-analysis does not include haploid abundance data of this species, we suspect, following our findings on the haploid–diploid paired species, that the haploid phase of E. huxleyi is also ecologically relevant. Previous studies suggest that the haploid life cycle phase of E. huxleyi has a different niche due to its streamlined metabolism (Rokitta et al., 2011) and variations in its response to bacterial (Mayers et al., 2016; Bramucci et al., 2018) and viral pressures (Frada et al., 2008). However, it should be noted that in some instances morphology rather than ploidy level seems to be the primary driver for observed differences in E. huxleyi (Frada et al., 2017). Overall, observations in the haploid stage of E. huxleyi are extremely limited due to difficulty of identifying the haploid phase with regular light microscopy, highlighting the need for developing new techniques to account for this potentially important life cycle stage. Further development of FISH (Campbell et al., 1994) and COD-FISH (Frada et al., 2012) methodologies would be particularly relevant in this context.

4.2 Concluding remarks

Our data compilation provides insight into the distribution of heterococcolithophores and holococcolithophores but also highlights many gaps in the data distribution and our knowledge on coccolithophore ecology.

There is for instance a lack of SEM observations in the Pacific Ocean (two studies in this compilation). However, this is in part because existing data from the 1980s were not retrievable and because the low abundance of coccolithophores in this region means that it has been of low priority for time-intensive and costly re-sampling. In addition, there is a limited number of available SEM time series, which are particularly valuable due to the seasonal nature of these organisms and the importance of time in structuring coccolithophore niche. The patchiness of the data combined with the patchiness of coccolithophore blooms is a challenge for fully assessing marine ecosystem functioning and in providing global abundance estimations.

Nonetheless, from our compilation it is clear that holococcolithophores constitute a minor component of total coccolithophore abundance. This could be in part due to sampling bias, specifically temporal bias towards periods of high heterococcolithophore abundance. However, other factors such as the strong dominance of E. huxleyi, which has a naked haploid phase, and the limited biomass low-nutrient regions are able to sustain might also exert a significant influence. The low contribution of holococcolithophores is interesting and raises the question of which physiological traits make heterococcolithophores generally more successful in the modern ocean.

Aside from limitations of in situ measurements, size and POC and PIC measurements of paired heterococcolithophore and holococcolithophore species are sparse, particularly for holococcolithophores.

Such measurements are needed for global organic carbon and carbonate production estimates, which are critical for biogeochemical estimates, including Earth system model studies. Models that could then be used to contextualize in situ observations in a biogeochemical context and test responses to environmental pressures presented by anthropogenic CO₂ emissions. Modelling approaches could furthermore be used to investigate drivers of distribution trends are that difficult to acquire with in situ measurements, such as the role of competition with other phytoplankton and the influence of top-down control on distribution trends, both of which have been shown to be important drivers of coccolithophore distribution in previous studies (Monteiro et al., 2016; Nissen et al., 2018).

A pertinent environmental driver not covered in our meta-analysis due to limited data is the influence of carbonate chemistry within the haploid–diploid niche. As the haploid and diploid phases of coccolithophores vary in their calcification status, they may therefore show different responses to carbonate chemistry. A study by Triantaphyllou et al. (2018) for instance found that holococcolithophores increased abundance in low pH waters. If this holds true on a global level, and holococcolithophores inhabit lower pH waters in terms of their niche, this would have important implications in the context of ocean acidification. In particular because meta-analysis (Ridgwell et al., 2009; Krumhardt et al., 2017) and modelling (Ridgwell et al., 2007; Krumhardt et al., 2019) suggest a shift towards lower global calcification rates in response to ocean acidification and warming. It should, however, be noted that the response of heterococcolithophores to ocean acidification is both strain and species dependent (Langer et al., 2006, 2009; Meyer and Riebesell, 2015), and global calcification rates might be more impacted by shifts in species composition rather than individual responses (Ridgwell et al., 2009). Furthermore, contradicting evidence suggesting increased coccolithophore abundance in response to higher CO₂ has been noted in situ (Rivero-Calle et al., 2015)

Finally, additional experiments on the numerical response of heterococcolithophores and holococcolithophores to various environmental drivers such as those performed on E. huxleyi would allow a better understanding of individual environmental pressures and will furthermore be highly valuable for future modelling approaches. In this context a better understanding of the triggers of phase transition would additionally be highly desirable, as the lack of haploid–diploid pairs of the same strain limits genomic approaches.

5 Conclusions

Our analysis shows that holococcolithophores constitute a minor proportion of total coccolithophore abundance (≈ 2 %–15 %) and constitute about ≈ 5 %–30 % of total paired coccolithophore abundance depending on the season.

Furthermore, our study shows that heterococcolithophores and holococcolithophores have contrasting environmental preferences and that therefore the haplo-diplontic life cycle expands the niche volume coccolithophores can inhabit by ≈17 %. Although our findings are limited to holococcolith-forming species, lab studies suggest similar patterns are likely to be observed for other coccolithophore species such as E. huxleyi, and this raises the question of how much the haploid phase of this species contributes to global coccolithophore abundance.

These results highlight the need to include haploid cells into coccolithophore studies in the context of environmental studies, modelling approaches, and physiological studies. We limit our understanding of these organisms by only focusing on one life cycle phase, particularly in the context of coccolithophore response to climate change, as increased stratification in a warming climate may favour the haploid life cycle of coccolithophores.