Late Neogene evolution of modern deep-dwelling plankton

The fossil record of marine microplankton provides insights into the evolutionary drivers which led to the origin of modern deep-water plankton, one of the largest components of ocean biomass. We use global abundance and biogeographic data, combined with depth habitat reconstructions, to determine the environmental mechanisms behind speciation in two groups of pelagic microfossils over the past 15 Myr. We compare our microfossil datasets with water column profiles simulated in an Earth system model. We show that deep-living planktonic foraminiferal (zooplankton) and calcareous nannofossil (mixotroph phytoplankton) species were virtually absent globally during the peak of the middle Miocene warmth. The evolution of deep-dwelling planktonic foraminifera started from subpolar–mid-latitude species, during late Miocene cooling, via allopatry. Deep-dwelling species subsequently spread towards lower latitudes and further diversified via depth sympatry, establishing modern communities stratified hundreds of metres down the water column. Similarly, sub-euphotic zone specialist calcareous nannofossils become a major component of tropical and sub-tropical assemblages during the latest Miocene to early Pliocene. Our model simulations suggest that increased organic matter and oxygen availability for planktonic foraminifera, and increased nutrients and light penetration for nannoplankton, favoured the evolution of new deep-water niches. These conditions resulted from global cooling and the associated increase in the efficiency of the biological pump over the last 15 Myr.

1 Introduction

The biodiversity of planktonic and nektonic organisms is difficult to explain, given the uniform character and vastness of pelagic environments where genetic isolation seems difficult to maintain (Norris, 2000). Planktonic microorganisms with mineralized shells have often been used as a model to study the mode and tempo of species origination in the open ocean, due to the abundance, widespread distribution, and temporal continuity of their fossil record (e.g. Pearson et al., 1997; Norris, 2000; Bown et al., 2004; Ezard et al., 2011; Norris et al., 2013). Because of the great fossilization potential of their calcium carbonate tests across much of the global ocean, their relatively simple and well-established taxonomy, and highly resolved biostratigraphy, planktonic foraminifera and calcareous nannofossils are amongst the most thoroughly studied. Planktonic foraminifera are heterotrophic zooplankton, with different species specialized to feed on different types of food, ranging from other plankton to sinking detritus. In the modern ocean, planktonic foraminifera live stratified across a range of depths spanning from the surface to hundreds of metres down the water column (Rebotim et al., 2017; Meilland et al., 2019). Properties such as food quantity and quality, oxygen, light, and pressure all change markedly across the first few hundreds of metres of the ocean. Depending on such down-column variability in environmental conditions, planktonic foraminifera can actively control their living depth of preference, which remains relatively stable during their adult life stage (Hull et al., 2011; Weiner et al., 2012; Rebotim et al., 2017; Meilland et al., 2019). A key advantage of using planktonic foraminifera for evolutionary studies is the ability to extract ecological information from their shell chemistry. This provides invaluable information about species–specific functional ecology (e.g. feeding strategy) and habitat preferences (e.g. surface versus deep waters), which, in combination with biogeographic, taxonomic, biometric, and stratigraphic data, have often been used to infer speciation and extinction mechanisms (Norris et al., 1993, 1994; Pearson et al., 1997; Hull and Norris, 2009; Pearson and Coxall, 2012; Woodhouse et al., 2021) and reconstruct phylogenetic relationships (Aze et al., 2011).

Calcareous nannoplankton also have a highly resolved and continuous fossil record; they are the most abundant microfossils in oceanic pelagic sediments, and similar to planktonic foraminifera, their spatial distribution ranges from tropical to subpolar latitudes (Poulton et al., 2017). In the modern ocean, they also live stratified in the water column, with species adapted to euphotic waters and species adapted to living deeper (Poulton et al., 2017). In contrast to planktonic foraminifera, nannoplankton are predominantly autotrophic, performing photosynthesis in water where light penetration is sufficient, although there is evidence for heterotrophy (mixotrophic behaviour) in some extant (Godrijan et al., 2020) and fossil (Gibbs et al., 2020) species. In euphotic waters, organic matter production from nannoplankton is at the base of pelagic food chains and of the functioning of the ocean biological carbon pump. Taxonomic, biometric, and stratigraphic data have been used to establish phylogenetic relationships between fossil nannoplankton species (Young and Bown, 1997).

Little emphasis has been given to the long-term drivers of evolutionary patterns observed in fossil plankton from species to phylum level, although, more recently, a broad connection with changing climate and ocean properties has been suggested (e.g. Ezard et al., 2011; Norris et al., 2013; Frass et al., 2015; Henderiks et al., 2020; Lowery et al., 2020). Boscolo-Galazzo et al. (2021) showed that, over the last 15 Myr, the remineralization of particulate organic carbon (POC) in surface waters declined markedly, driven by climate and ocean cooling (Kennett and Von der Borch, 1986; Kennett and Exon, 2004; Cramer et al., 2011; Zhang et al., 2014; Herbert et al., 2016; Sosdian et al., 2018; Super et al., 2020), increasing the efficiency of the ocean biological pump in delivering organic matter at depth. Such a mechanism was key in promoting the evolution of life in deep waters, allowing the development of the modern twilight zone ecosystem (Boscolo-Galazzo et al., 2021). The goal of this study is to combine the fossil record of two ecologically complementary calcareous microplankton groups seldom analysed together, planktonic foraminifera and nannoplankton, and, together with model simulations, help disentangle the evolutionary drivers of modern deep-dwelling plankton. We use the planktonic foraminiferal dataset from Boscolo-Galazzo et al. (2021) and extend our analysis to calcareous nannofossils in coeval sediment samples to assess their abundance and distribution pattern. We compare the results from these two groups and contrast them against time and site-specific model water column profiles for POC and oxygen (O₂) obtained from the cGENIE Earth system model. Furthermore, using stable isotopes, depth habitat reconstructions, abundance, and biogeography data, we reconstruct the speciation mechanisms which led to the evolution of modern deep-dwelling planktonic foraminiferal species.

2 Methods

2.1 Planktonic foraminifera

In this study, we focus on the deep-dwelling groups of macroperforate planktonic foraminifera of the hirsutellids, globorotaliids, truncorotaliids, and globoconellids, which, in the modern ocean, calcify and live mostly in the twilight zone of the ocean, i.e. between 200–1000 m (Birch et al., 2013; Rebotim et al., 2017), and have a more complete fossil record than deep-dwelling microperforate planktonic foraminifera (Kennett and Srinivasan, 1983).

Figure 1 Depth habitat reconstructions for middle Miocene to present planktonic foraminiferal species at the investigated sites. Surface dwellers (species living at depths shallower than 200 m) are indicated with a grey dot, and deeper species are indicated with coloured symbols. The relative size of the symbols represents the size fractions of the sample. Reproduced from Boscolo-Galazzo et al. (2021).

Planktonic foraminiferal data and depth habitat reconstructions (Fig. 1) are from Boscolo-Galazzo et al. (2021). They were obtained from the globally and latitudinally distributed DSDP (Deep Sea Drilling Project), ODP (Ocean Drilling Program), and IODP (Integrated Ocean Drilling Program and International Ocean Discovery Program) sites and from cores drilled onshore and offshore Tanzania, which are all characterized by abundant calcareous microfossils (Boscolo-Galazzo et al., 2021). The work focused on seven target ages (15, 12.5, 10, 7.5, 4.5, 2.5, and 0 Ma in the Holocene). To avoid sample aliasing, bulk sediment stable isotopes were measured on an average of 10 samples per target age at each site. The sample displaying mean oxygen stable isotope values was chosen for subsequent analyses (Boscolo-Galazzo et al., 2021). The taxonomy follows the work of Berggren (1977), Kennett and Srinivasan (1983), Bolli et al. (1989), Scott et al. (1990), Berggren (1992), Pearson (1995), Majewski (2010), Fox and Wade (2013), Spezzaferri et al. (2015), Wade et al. (2018), and Lam and Leckie (2020a), with phylogenetic genus names from Aze et al. (2011).

Ages were determined based on biostratigraphic analysis, mostly following the biozonation scheme by Wade et al. (2011).

Foraminiferal picking for stable isotope measurements were conducted from the following size fractions: 180–250, 250–300, and 300–355 µm (Boscolo-Galazzo et al., 2021). Stable isotopes were measured on an average of 15 different species per sample, using ∼ 25 specimens for common species, and as many specimens as possible for rare species. Stable isotopes were measured at Cardiff University. Stable isotope results are shown in Figs. S1 to S9 in the Supplement. Only data from the largest of the three measured size fractions are shown when data from more than one size fraction are available. Data from size fractions other than those above are shown only when a species did not occur within the preferred size interval. Foraminiferal abundance counts were carried out in two size fractions, 180–250 and > 250 µm, counting up to 300 specimens in each. Total abundances were derived by summing up the abundances from these two size fractions.

Boscolo-Galazzo et al. (2021) reconstructed the planktonic foraminiferal depth habitat (Fig. 1) using a combined model–data approach, solving the palaeotemperature equation of Kim and O'Neill (1997) for each data point, using measured foraminiferal δ¹⁸O, global ice volume estimates, and the cGENIE modelled salinity field to determine local water δ¹⁸O and then using the model temperature–depth curve to determine depth. The full method is described in Boscolo-Galazzo et al. (2021).

2.2 Calcareous nannofossils

Quantitative calcareous nannofossil data were collected from the same samples as used for planktonic foraminiferal analysis or, when this was not possible, stratigraphically adjacent samples (Table S1 in the Supplement). A cascading count technique was used to maximize the nannofossil diversity recovery and the quantification of low abundance species (Styzen, 1997). Nannofossils were counted per field of view (FOV) until a minimum of 400 specimens were achieved for each sample. However, if a high-abundance species exceeded an average of 25 specimens per FOV, it was excluded from subsequent counts in that sample, and its abundance was scaled up based on its average abundance and the total number of FOVs counted. Only specimens directly counted contributed to the minimum count threshold of 400 specimens. An additional scan of two slide transects were undertaken to record rare species not observed during the extended count and are included in the total species richness and diversity analyses. Samples for nannofossil analysis were prepared using the smear slide technique (Bown and Young, 1998). Calcareous nannofossils were observed using both plane-polarized light (PPL) and cross-polarized light (XPL) on a ZEISS Axioscope light microscope at × 1000 magnification. The identification and taxonomy used herein follow Young et al. (1997) and are coherent with the recent Neogene calcareous nannofossil taxonomy (Ciummelli et al., 2016; Bergen et al., 2017; Blair et al., 2017; Boesiger et al., 2017; Browning et al., 2017; de Kaenel et al., 2017).

2.3 Plankton ecogroups

In order to compare the datasets obtained from the planktonic foraminiferal and nannofossil analysis, we grouped the species into ecogroups based on depth habitat preferences. Planktonic foraminiferal ecogroups are defined based on palaeodepth habitat reconstructions from Boscolo-Galazzo et al. (2021). The euphotic zone ecogroups includes species with an average depth habitat shallower than 200 m (the bottom of the euphotic zone), and the twilight zone ecogroup includes species with an average depth habitat coinciding with the twilight zone (200–1000 m). The twilight zone ecogroup is largely composed of species within the globoconellids, the Globorotalia merotumida–tumida lineage, the hirsutellids, and the truncorotaliids but also includes species from other genera, such as Globigerinella calida, Globorotaloides hexagonus, G. variabilis, and Pulleniatina obliquiloculata. Dentoglobigerina venezuelana has a changeable depth habitat through time (Matsui et al., 2016; Wade et al., 2018); following the depth habitat reconstructions from Boscolo-Galazzo et al., (2021) it was grouped as euphotic zone species for target ages 15, 12.5, 7.5, and 4.5 Ma and as twilight zone species for target age 10 Ma. Species were excluded from the grouping when they are known to have a marked seasonality in abundance and depth habitat (Globigerina bulloides, Globigerinella praesiphonifera and the neogloboquadrinids; e.g. Jonkers and Kucera, 2015; Greco et al., 2019), and if they were too rare, and depth habitat reconstruction was not possible (Candeina nitida).

A total of three ecogroups for calcareous nannofossils are used, namely upper euphotic, lower euphotic, and sub-euphotic. The upper-euphotic group is represented by Discoaster spp., Rhabdosphaera xiphos, Reticulofenestra spp., and Gephyrocapsa spp. (excluding G. ericsonii); the lower-euphotic ecogroup contains Rhabdosphaera clavigera, Gephyrocapsa ericsonii, and Ceratolithus spp. Finally, the sub-euphotic ecogroup includes Florisphaera profunda and Calciosolena murrayi. Because species-specific stable isotope measurements and depth habitat reconstructions are difficult for calcareous nannofossils, species depth habitat preference was assigned based on the literature (Poulton et al., 2017; Tangunan et al., 2018). In particular, Poulton et al. (2017) described vertically separated coccolithophores communities sampled during a meridional cruise in the Atlantic Ocean. Here we use their criteria for assigning nannofossil species into ecogroups, whereby, in the upper-euphotic zone ecogroup, we include species found to live in waters with > 10 % surface irradiance. In the lower-euphotic zone ecogroup, we group species found to live in waters with 10 %–1 % irradiance, and in the sub-euphotic zone ecogroup, we group species found to live in waters with < 1 % irradiance, i.e. too low to support photosynthesis (Poulton et al., 2017). Discoaster become extinct in the early Pleistocene; therefore, its depth habitat remains under debate, as the group has no extant relative (Schueth and Bralower, 2015; Tangunan et al., 2018). However, geochemical evidence from oxygen isotope values of Discoaster and planktonic foraminifera (Globorotalia menardii, Dentoglobigerina altispira, and Globigerinoides obliquus) reveal comparable values and signify that Discoaster likely inhabited the upper-euphotic zone (Minoletti et al., 2001).

Figure 2 Foraminiferal and nannofossil ecogroup abundance at Sites 242 and U1482.

Figure 3 Foraminiferal and nannofossil ecogroup abundance at Sites 871/872 and U1490.

Figure 4 Foraminiferal and nannofossil ecogroup abundance at Sites U1338 and U1489.

Figure 5 Foraminiferal and nannofossil ecogroup abundance at Sites 516 and 1138.

Figure 6 Relative abundance and diversity indexes of planktonic foraminiferal and nannofossil ecogroups.

For each target age, the relative abundance of ecogroups was calculated by summing up the abundance counts of all the individual species pertaining to an ecogroup at each site; hence, ecogroup abundance data represent global mean values. For both nannofossils and planktonic foraminifera, the percentage of each ecogroup per time bin was converted into pie charts (Figs. 2–5). Diversity indexes for both the foraminiferal and nannofossil ecogroup were calculated using the statistical software Past (Hammer et al., 2001; Fig. 6).

2.4 cGENIE model

We extracted the model output for the particulate organic carbon (POC) and oxygen concentration from the cGENIE simulations for each of the seven target ages (Zenodo: https://doi.org/10.5281/zenodo.4469673 and https://doi.org/10.5281/zenodo.4469678, Ridgwell et al., 2021a, b), as described fully in Boscolo-Galazzo et al. (2021). To facilitate a general discussion of near-surface changes, we divided the data latitudinally by calculating the arithmetic mean for low latitudes (< 16^∘ latitude), two mid-latitude bands (mid 1 at 16 to 40^∘; mid 2 at 40 to 56^∘) and high latitudes (> 56^∘). The cGENIE simulations take changing boundary conditions into account, including CO₂ forcing, bathymetry, and ocean circulation (Crichton et al., 2021a). The model's ocean biological carbon pump is temperature dependent, where temperature affects both nutrient uptake rates at the surface and remineralization rates of sinking particulate organic matter down the water column (Crichton et al., 2021b).

3 Results

3.1 Plankton ecogroups

For both calcareous nannoplankton and planktonic foraminifera, the variation in the abundance of euphotic zone and deeper-dwelling ecogroups show global patterns recognized across sites. Additionally, both group indicate a long-term directionality towards increased abundance of deep-dwelling ecotypes. Among planktonic foraminifera, the twilight zone ecogroup increases in abundance through time, starting at 7.5 Ma (Fig. 6). The relative abundance of the twilight zone ecogroup in the middle Miocene is 15 %, and it increases to ∼ 30 % in the Holocene time slice (Fig. 6). The average abundance of the euphotic zone species ecogroup in the middle Miocene is 85 %, and it decreases through time until it reaches 60 % in the Holocene (Fig. 6). In the twilight zone ecogroup, we observe an increase in the total number of species from about 1–2 species at 15 Ma to 14 species in the Holocene (Fig. 6). In the middle Miocene, this group comprised one-sixth of the total number of species in our samples, while in the Holocene it represents almost the half. All the diversity indexes show a late Miocene to Holocene increasing trend for the twilight zone ecogroup (Fig. 6).

Calcareous nannofossil assemblages are dominated by the upper-euphotic ecogroup from 15 to 10 Myr at all sites (Figs. 2–5). At 7.5 Ma, the sub-euphotic ecogroup first becomes a significant component of assemblages at Indian Ocean sub-tropical Site U1482 and, to some extent Site 242, but it is not until the 4.5 Ma time slice that the sub-euphotic ecogroup becomes a significant component of assemblages at the majority of locations (Sites 516, 871/872, 242, U1338, U1482, and U1489; Figs. 2–5). By the Holocene time slice, coccoliths of sub-euphotic species are dominant at most locations, except at eastern equatorial Pacific Site U1338 (Figs. 6 and 4). At the southern high latitude Site 1138, there is no significant contribution from coccoliths of either lower-euphotic or sub-euphotic species at any point, although there are no data from the Pliocene to Holocene time slices at this location (Fig. 5). Global average compositions of calcareous nannofossil assemblages reflect the changes noted above, with a marked and rapid decline in the relative contribution of the upper-euphotic ecogroup and a corresponding increase in the sub-euphotic zone ecogroup through the Pliocene and to Holocene (Fig. 6).

Figure 7 Abundance and biogeography of planktonic foraminiferal species at 15, 12.5, and 10 Ma. In the pie charts, the twilight zone planktonic foraminiferal species are in bold. Sites where twilight zone planktonic foraminifera were found are highlighted in green in the maps. The crossed square symbols in the maps indicate that the time interval of interest was not recovered for a given site. Continental configuration follows the Ocean Drilling Stratigraphic Network (GEOMAR; Kiel, Germany).

Figure 8 Abundance and biogeography of planktonic foraminiferal species at 7.5 Ma. In the pie charts, twilight zone planktonic foraminiferal species are in bold. Sites where twilight zone planktonic foraminifera were found are highlighted in green on the maps. The crossed square symbols in the maps indicate that the time interval of interest was not recovered for a given site. Continental configuration follows the Ocean Drilling Stratigraphic Network (GEOMAR; Kiel, Germany).

3.2 Planktonic foraminiferal deep-dwelling species: depth habitat, abundance, and biogeography

3.2.1 Hirsutellids

The only hirsutellid species occurring in our Miocene samples is Hirsutella scitula. At 15 Ma, this morphospecies is common only at Site 1138 (∼ 8 %), occurs sporadically at Site 516 (< 1 %), and is absent at the other investigated sites (Fig. 7). Oxygen isotopes range from 0.5 ‰ to 1.3 ‰ (Figs. S1 and S2). Depth habitat reconstructions for Site 1138 are unattainable from δ¹⁸O data due to the overprinting effect of subpolar front shifts at this location, but a habitat reconstruction at Site 516 suggests a palaeodepth habitat shallower than 200 m. By 12.5 Ma, H. scitula appears at Sites U1338 and U1489 in very low abundance (< 0.5 %; Fig. 7). At Site U1489, the species was so rare that it was encountered when picking for stable isotopes and no more when counting for species abundances, despite the use of different splits of the residue. No differences were observed between 12.5 and 10 Ma in the biogeography of H. scitula (Fig. 7). However, by 7.5 Ma, H. scitula occurs at all our low-latitude sites (Fig. 8), with oxygen isotopes ranging from −0.5 ‰ to 2.0 ‰ (Figs. S3–S9), which, according to depth habitat reconstructions, translates to 250 and 500 m water depth (Boscolo-Galazzo et al., 2021). This is similar to that of Globorotaloides hexagonus (Fig. 1), which is the only twilight zone dweller we observed at tropical sites at 15 Ma, displaying stable isotopes ranging from 0 ‰ to 1 ‰, which translates to depths around 300–500 m. In the late middle Miocene, the stable isotope values of G. hexagonus increase to 2 ‰–2.5 ‰ (Figs. S3–S9). Similarly, the oxygen isotope values of tropical H. scitula increased through time, reaching 2 ‰–3 ‰ in the youngest target ages. In line with this, the reconstructed depth habitat of H. scitula and G. hexagonus increases through time in a stepwise fashion, and in the Holocene it reaches down to 800–1500 m (Fig. 1; Boscolo-Galazzo et al., 2021). Hirsutella scitula becomes gradually more common at low-latitude sites through the Miocene–Pliocene, although it never becomes abundant. In our record, Hirsutella margaritae and H. theyeri first appear in the early Pliocene at a depth between 200–300 m (Fig. 1; oxygen isotopes range from −1 ‰ to −0.5 ‰), similar to that of H. hirsuta (oxygen isotopes range from 0 ‰ to 1 ‰) when it first appears in the Holocene (Figs. S3–S9).

Figure 9 Abundance and biogeography of planktonic foraminiferal species at 4.5 Ma. In the pie charts, twilight zone planktonic foraminiferal species are in bold. Sites where twilight zone planktonic foraminifera were found are highlighted in green in the maps. The crossed square symbol in the maps indicate that the time interval of interest was not recovered for a given site. Continental configuration follows the Ocean Drilling Stratigraphic Network (GEOMAR; Kiel, Germany).

Figure 10 Abundance and biogeography of planktonic foraminiferal species at 2.5 Ma. In the pie charts, twilight zone planktonic foraminiferal species are in bold. Sites where twilight zone planktonic foraminifera were found are highlighted in green in the maps. The crossed square symbols in the maps indicate that the time interval of interest was not recovered for a given site. Continental configuration follows the Ocean Drilling Stratigraphic Network (GEOMAR; Kiel, Germany).

3.2.2 Truncorotaliids

The earliest appearances of Truncorotalia crassaformis in our record corresponds to our 4.5 Ma time slice at Site 1138 in the Indian Ocean sector of the Southern Ocean, where it represents > 20 % of the assemblage, and in coeval sediments at mid-latitude Site 516 in the southwest Atlantic, where it represents ∼ 9 % of the assemblage (Fig. 9), with oxygen isotopes ranging from 1.5 ‰ to 3.0 ‰ (Figs. S1 and S2). At Site 516, we observe the co-occurrence of T. oceanica (∼ 14 %) and T. crassaformis in the early Pliocene (4.5 Ma) and of T. viola (∼ 5 %) and T. crassaformis (∼ 13 %) in the late Pliocene (2.5 Ma; Fig. 10). We did not observe T. oceanica and T. viola anywhere else. T. oceanica and T. crassaformis display almost overlapping δ¹⁸O values and depth habitat (Fig. S2), but T. crassaformis has 0.5 ‰ lower δ¹³C values. Oxygen stable isotope data (1 ‰; Fig. S2) and habitat reconstructions for Site 516 indicate that a subsurface habitat (> 200 m) was already occupied by T. crassaformis at the beginning of its evolutionary history (Fig. 1). The late Pliocene appearance of T. viola at Site 516, which differs from T. crassaformis in having a more convex umbilical side, a triangular outline, and a subacute profile, is associated with a shift to more positive oxygen isotope values of T. crassaformis (1.5 ‰) and to slightly greater depths (Fig. 1).

Figure 11 Abundance and biogeography of twilight zone planktonic foraminiferal species at 0 Ma. In the pie charts, twilight zone planktonic foraminiferal species are in bold. Sites where twilight zone planktonic foraminifera were found are highlighted in green in the maps. The crossed square symbols in the maps indicate that the time interval of interest was not recovered for a given site. Continental configuration follows the Ocean Drilling Stratigraphic Network (GEOMAR; Kiel, Germany).

We find T. crassaformis by the late Pliocene at our investigated tropical and sub-tropical sites (2.5 Ma; Sites U1338, U1489, 872, U1490, 242, and U1482; Figs. 10 and 11), with oxygen isotope values ranging from 1.0 ‰ to 2.0 ‰, which translate to depth habitats of 400–600 m (Boscolo-Galazzo et al., 2021; Fig.1). The appearance of T. crassaformis at our low-latitude sites is coeval with the appearance in our record of Truncorotalia tosaensis, morphologically transitional between T. crassaformis and T. truncatulinoides (Lazarus et al., 1995; Fig. 10). Truncorotalia tosaensis displays oxygen isotopes values ranging from 0 ‰ to 0.5 ‰ (Figs. S1–S9), which translate to 300–350 m depth (Fig. 1). Consistent with earlier findings (Jenkins and Srinivasan, 1986; Lam and Leckie, 2020b; Lazarus et al., 1995), we record the first occurrence of T. truncatulinoides in the late Pliocene in the southwestern Pacific (2.5 Ma; Site U1482), and only later in the North Pacific (0 Ma; Site 872), Indian Ocean (0 Ma; Site 242), and South Atlantic (0 Ma; Site 516; Fig. 11). Truncorotalia truncatulinoides records oxygen isotope values ranging from −1 ‰ to 2 ‰, which are more negative than coeval T. crassaformis (Figs. S3–S9). Truncorotalia truncatulinoides, although reported in the modern tropical ocean as one of the species living at the greatest depths, occupies a shallower depth habitat than T. crassaformis when it first appears in our tropical to sub-tropical records (2.5 Ma).

3.2.3 Globorotaliids

With globorotaliids, here we refer to the Globorotalia merotumida–tumida lineage composed by G. merotumida, G. plesiotumida, G. tumida, and G. ungulata (Kennett and Srinivasan, 1983). This group first appears with Globorotalia plesiotumida in our 10 Ma time slice at Site 871 (Fig. 7). At all the investigated low-latitude sites, we find G. tumida by 4.5 Ma, with abundances between 2 % and 24 % (Fig. 9). Globorotalia plesiotumida co-occurs with G. tumida only at Sites 872 and 242, corresponding to our 4.5 Ma time slice (Figs. 9 and S9). In our records, G. tumida consistently displays oxygen isotope values between −1 ‰ and 0 ‰ (Figs. S3–S9) and an average depth habitat around 250 m, with a shallowest occurrence at 50 m and a deepest occurrence around 600 m (Fig. 1). Similar oxygen isotope values and depth habitat preference are recorded for G. plesiotumida and G. ungulata (Figs. 1 and S3–S9).

3.2.4 Globoconellids

In our records, Globoconella miozea is a dominant component of planktonic foraminiferal assemblage at Site 1138 at 15 Ma (65 %) and occurs in moderate abundance at Site 516 (5 %; Fig. 7). The distribution of globoconellids appears restricted to southern high to mid-latitudes during the middle Miocene and the late Miocene to Pliocene. At Site 1138, globoconellids decrease in abundance through time, with Globoconella panda being the only late Miocene (7.5 Ma) species (< 1 %; with δ¹⁸O of 3 ‰), followed only by Globoconella inflata in the Holocene (4.2 %; with δ¹⁸O of 3.5 ‰; Fig. 11). On the contrary, at Site 516, globoconellids increase in abundance through time, becoming a characteristic feature of the planktonic foraminifera assemblage, as noted in previous studies of this area (Berggren, 1977; Norris et al., 1994; Figs. 7–11). In the Holocene, G. inflata is most abundant at the mid-latitude Site 516 (22.7 %; with δ¹⁸O of 1 ‰) but also occurs in the sub-tropical (< 0.5 %; with δ¹⁸O of 0.9 ‰) and subpolar Indian Ocean (4.2 %; with δ¹⁸O of 3.5 ‰; Figs. 11 and S1–S3). Depth habitat reconstructions for the globoconellids show a deepening trend through time, although it is less marked compared to those of other deep-dwelling groups considered in this study (Boscolo-Galazzo et al., 2021; Fig. 1). At Site 516, the depth habitat is just above 200 m for middle Miocene Globoconella miozea (δ¹⁸O 0.9 ‰), just below 200 m for late Miocene G. miotumida and G. conomiozea (δ¹⁸O 1 ‰), and around 350 m for late Pliocene Globoconella puncticulata (δ¹⁸O 1.6 ‰), the precursor of G. inflata, which shows a similar depth habitat at this site (δ¹⁸O 1 ‰–1.2 ‰; Figs. S2 and 1). At Site 242, the average depth reconstruction for the Holocene G. inflata is of ∼ 450 m (δ¹⁸O 0.9 ‰; Fig. S3), which is similar to that of T. crassaformis and T. truncatulinoides (Boscolo-Galazzo et al., 2021; Fig. 1).

4 Discussion

4.1 Evolution of a deep-plankton ecological niche linked to late Neogene cooling

We observe a trend of increasing ecological importance of deep-dwelling species in both calcareous nannofossil and foraminiferal communities from the late Miocene to the present. The foraminiferal twilight zone ecogroup shows a marked increase in abundance and diversity at 7.5 Ma, at the same time as the first significant appearance of the sub-euphotic ecogroup within coccolithophore assemblages and coinciding with a possible acceleration of global cooling (Cramer et al., 2011; Crichton et al., 2021a). Fossil deep-dwelling coccoliths are dominated by one species, Florisphaera profunda, which is often very abundant in Plio–Pleistocene sediments. This makes our sub-euphotic ecogroup a low diversity–high abundance assemblage, whose origin significantly impacts the relative abundance balance between ecogroups but with little change in diversity metrics towards the modern age. Although there are morphological variants of F. profunda in the modern oceans, potentially representing sub-species or pseudo-cryptic species (Quinn et al., 2005), these are not distinguished in fossil assemblages, and documenting their divergence times requires either further molecular genetic or detailed morphological analyses. The one clear signal in our diversity analyses is a late Miocene–early Pliocene peak in upper-euphotic species richness, followed by a marked decline through the late Pliocene to the Holocene. The late Miocene–Pliocene peak diversity is present in previous global compilations of total nannofossil diversity (Bown et al., 2004; Lowery et al., 2020), but here, we show that this signal is driven first by a diversification and then by progressive extinction almost entirely within the upper-euphotic taxa.

Modern planktonic foraminifera evolved in two main diversification pulses in the middle Miocene (16–14 Ma) and during the late Miocene–Pliocene transition (6–4.5 Ma; Kennett and Srinivasan, 1983; Kucera and Schönfeld, 2007). Our species abundance and diversity data show that this diversification was mostly driven by the origin of lineages of deep/subsurface dwelling species (Figs. 7–11 and 6). Diversity among euphotic zone species remained constant from the middle Miocene to the early Pliocene and then declined (Fig. 6). This pattern, similar to calcareous nannofossils, may explain early records pointing to a decrease in Pliocene to the present planktonic foraminiferal diversity (Wei and Kennett, 1986).

The observed evolutionary patterns in planktonic foraminifera and coccolithophores can be explained by the development of environmental conditions favourable to deep-living plankton with cooling. For the deep-dwelling planktonic foraminifera, survival requires the availability of food at depth and, with the exception of few species adapted to oxygen minimum zones (Davis et al., 2021), the absence of severe oxygen depletion. In the published literature, there is a general tendency for planktonic foraminiferal δ¹⁸O values to be tightly grouped during times of warm climate, becoming more spread out during cooling episodes such as, for instance, the transition from the mid to late Cretaceous (Ando et al., 2010) and early to middle Eocene (John et al., 2013). In light of our findings, such tight δ¹⁸O values may indicate times when food and oxygen availability at depth were limited, allowing planktonic foraminifera to live at shallower depths only. Because our data are restricted to certain time slices and locations, we only capture a part of the overall pattern for the Neogene, but other examples of depth-related evolution in the Neogene include the Fohsella peripheroronda–F. fohsi lineage in the middle Miocene (Hodell and Vayavananda, 1993; Norris et al., 1993) and the appearance of various deep-dwelling digitate species in the Plio–Pleistocene (Coxall et al., 2007).

Deep-dwelling coccolithophores, most notably Florisphaera profunda, require dissolved macronutrients (N and P) and at least some degree of light penetration (Quinn et al., 2005; Poulton et al., 2017). The requirement for light penetration to depth is most clearly shown by the absence of F. profunda beneath high surface productivity regions, where nutrients are not limited at depth but light is rapidly attenuated by more abundant mixed layer microplankton (Beaufort et al., 1997). For coccolithophores, the cooling-driven shift from nutrient recycling within to below the mixed layer may have provided the ecological driver for species to live in deeper, more nutrient-rich waters but, as mixed layer waters cleared, also allowed the irradiance necessary for photosynthesis to penetrate to these new deeper habitats. Additionally, the capability of coccolithophores to absorb carbon and nutrients from seawater under low light conditions (Godrijan et al., 2020, 2022) may have also aided in the occupation of new deep-water niches.

Figure 12 cGENIE model output for changes in particulate organic carbon (POC) flux and dissolved oxygen in near-surface ocean waters from the middle Miocene to the present, with a temperature-dependent biological carbon pump. The inset map shows the modelled present surface ocean temperatures. Depths are the middle of cGENIE's top three ocean layers.

Our model outputs from cGENIE are consistent with this interpretation of organic matter export (POC at 40 m in Fig. 12) being reduced with cooling, suggesting an overall decrease in primary productivity at all four latitudinal bands considered in this study. Fewer particles in surface waters would have allowed greater light penetration (Fig. 12). At the same time, the model indicates enhanced organic matter delivery at > 200 m with cooling. Greater organic matter delivery below the euphotic zone suggests a deeper remineralization depth and increased dissolved nutrient availability at depth. This is most clearly shown in the modelled low- and mid-latitude near-surface waters. Oxygen availability also increased, particularly below 100 m depth in low latitudes and further down the water column below 200 m depth (Fig. 12).

Nonetheless, planktonic foraminifera and nannoplankton have distinct trophic statuses (zoo- versus phytoplankton). Furthermore, coccolithophore species require light for their dominantly photosynthetic mode of life. In our data, such differences between the requirements of zoo- and phytoplankton deep-dwellers is clearly observed in the biogeographic patterns. Sub-euphotic coccolithophores are consistently more abundant in low nutrient sub-tropical locations (e.g. DSDP Site 242 and IODP Site U1482; Fig. 2). The end-member of this biogeographic difference is ODP Site 1138 in the Southern Ocean. Here, twilight foraminifera dominate most time intervals, presumably due to high export production. However, lower-euphotic and sub-euphotic coccolithophores are effectively excluded due to turbid, high-productivity surface waters (Fig. 5). This pattern is also supported through a comparison of high to low productivity tropical sites in the Holocene; the eastern Equatorial Pacific Site U1338 has abundant twilight foraminifera but relatively low abundances of sub-euphotic coccolithophores (Fig. 4), whereas the more oligotrophic Site U1482 (Fig. 2) has lower abundances of twilight foraminifera and higher abundances of sub-euphotic coccolithophores.

Despite these ecological differences between zoo- and phytoplankton, we suggest that there is a shared environmental driver for the evolution of deep-dwelling coccolithophores and planktonic foraminifera linking the evolution of deep-dwelling specialists in each group. Efficient near-surface recycling of organic carbon in past warm climates, such as the middle Miocene (Fig. 12), precluded the occupation of the deep habitat for both groups by reducing both organic carbon transfer (food limitation for foraminifera and for foraminiferal prays such copepods) and light penetration (photosynthesis for coccolithophores) to depth (Figs. 2–6). Global cooling since the middle Miocene (Kennett and Von der Borch, 1986; Kennett and Exon, 2004; Cramer et al., 2011; Zhang et al., 2014; Herbert et al., 2016; Sosdian et al., 2018; Super et al., 2020) has, however, led to a decreased export of organic matter and a deepening of the mean organic matter remineralization depth, which in turn favoured the evolution of deep-water niches in planktonic foraminifera (Fig. 1) and nannoplankton via increased availability of organic matter, oxygen, and likely nutrients at depth (Fig. 12) and clearing of surface waters.

4.2 Mechanisms of speciation of deep-dwelling planktonic foraminifera

The hirsutellids gave rise to a large late Neogene radiation among planktonic foraminifera, leading to the origin of modern phyletic groups such as the menardellids, globoconellids, truncorotaliids, and the globorotaliids of the Globorotalia merotumida–tumida lineage (Kennett and Srinivasan, 1983; Scott et al., 1986; Aze et al., 2011). The majority of the modern representatives of these groups are lower-euphotic zone to twilight zone species. The hirsutellids originated about 18 Ma (Wade et al., 2011) from Globorotalia zealandica (Kennett and Srinivasan, 1983), with the first representative of the group being Hirsutella praescitula (Kennett and Srinivasan, 1983; Aze et al., 2011). Extant hirsutellids include H. scitula, H. hirsuta, and H. theyeri; genetic data available for H. hirsuta indicate a single genotype (Schiebel and Hembleben, 2017). Modern H. scitula and H. hirsuta are deep-water forms. Depth habitat reconstructions show H. scitula at the greatest water depth (Boscolo-Galazzo et al., 2021; Fig. 1), which is consistent with it being reported to have a deeper average depth habitat than H. hirsuta in the modern ocean (Birch et al., 2013; Stainbank et al., 2019), where it feeds on suspended organic matter (Itou et al., 2001). Hirsutella hirsuta has been reported to feed on dead diatoms and to predominantly live at depths around 250 m (Schiebel and Hemleben, 2017), consistent with our habitat reconstructions placing H. hirsuta and other hirsutellids shallower than H. scitula in the region between the bottom of the euphotic zone and the upper twilight zone. Our depth habitat reconstruction for H. scitula at Site 516 for the 15 Ma time slice indicates an initial euphotic zone depth habitat preference for this species. By the 7.5 Ma time slice, we find H. scitula at most of our low-latitude sites, at depths between 300 and 500 m (Fig. 1). We suggest that the spread of H. scitula from high–mid-latitudes towards the tropics after the middle Miocene warmth (Figs. 7–11) tracks the increasing availability of POC and oxygen at depth (Boscolo-Galazzo et al., 2021; Fig. 12), allowing this species to find food in deep tropical twilight zone waters and profiting from a new ecological niche. We suggest that moving from a surface to a deep-water habitat was, for H. scitula, an ecological innovation which allowed the species to move outside its high–mid-latitude areal, consistent with observations from the modern ocean documenting H. scitula dwelling at progressively deeper depth from high to tropical latitudes (Schiebel and Hemleben, 2017). The proceeding ocean cooling (and increasing efficiency of the biological pump) explains the stepwise depth habitat increase of H. scitula through time at tropical and sub-tropical sites (Fig.1) because more food became available at increasingly greater depth, in association with improved oxygen availability (Fig. 12), allowing the species to expand its habitat vertically rather than geographically. At low-latitude sites, this was accompanied by speciation, with the appearance of Hirsutella margaritae and H. theyeri in the early Pliocene and of H. hirsuta in the Holocene. These new species display depth habitats shallower than H. scitula at their appearance. Hence, we suggest that it was the occupation of a new deep-water habitat in tropical waters that triggered speciation from H. scitula through depth sympatry, i.e. genetic isolation attained in the same area but at different depths (Weiner et al., 2012). Earlier studies on speciation among planktonic foraminifera based on the fossil record highlighted a predominance of sympatric speciation, possibly linked to changes in ecology (Norris et al., 1993; Lazarus et al., 1995; Pearson et al., 1997). This has more recently been supported by genetic analysis, which reveals a consistent depth separation between intra-specific genotypes at a global scale, suggesting that depth sympatry could be a universal mechanism generating diversity among microplankton (Weiner et al., 2012).

Truncorotaliids start their evolutionary history at ∼ 4.31 Ma (Raffi et al., 2020), when Truncorotalia crassaformis splits from Hirsutella cibaoensis (Kennett and Srinivasan, 1983; Aze et al., 2011), which is one of the first new species originating from H. scitula after its spread to low latitudes (Kennett and Srinivasan, 1983). Our reconstructed geographic and temporal distribution for T. crassaformis suggests that the early Pliocene split from the hirsutellids happened in cold subpolar water. By the early Pliocene, Truncorotalia crassaformis was an abundant component of the subpolar assemblage at Site 1138 (> 17 %) and had already successfully spread to mid-latitudes (> 8 %) but does not occur in our low-latitude samples (Fig. 9). Although modern planktonic foraminifera have a high dispersal potential (Darling et al., 2000; Norris and De Vargas, 2000), our data suggest that the evolution of Truncorotalia crassaformis from the hirsutellids might have happened through allopatry, with the low-latitude population of H. cibaoensis becoming isolated from subpolar populations and eventually evolving into T. crassaformis (e.g. following the development of diachrony reproductive seasons in high- and low-latitude populations). Similar to H. scitula, from subpolar latitudes, Truncorotalia crassaformis appears to have subsequently spread to lower latitudes, occupying progressively deeper habitats (Figs. 1 and 7–11), and originating numerous daughter species, some of which are intermediate morphospecies with limited geographic distributions no longer extant today (Lazarus et al., 1995). At Site 516, we find Truncorotalia viola, with lighter δ¹⁸O and δ¹³C values than T. crassaformis (Fig. S2), pointing to a clear ecological differentiation. Together with the marked morphological differentiation between the two, this suggests that T. viola may have been a different biological species. Truncorotalia truncatulinoides, the most representative species of this group, appears to have originated from T. crassaformis at about 2.7 Ma in the tropical southwestern Pacific and subsequently spread in the global ocean (Dowsett, 1988; Lazarus et al., 1995; De Vargas et al., 2001). According to our reconstruction, depth sympatry associated with gradual morphological changes characterizes speciation among the truncorotaliids, as depth habitat deepening of the ancestor/precursor is clear in the transitions of T. crassaformis–T. viola and T. crassaformis–T. tosaensis–T. truncatulinoides in our record (Fig. 1). The possibility to colonize deeper water habitats may have led to progressive reproductive isolation between deeper and shallower populations of T. crassaformis, resulting in biological speciation. Depth sympatry as a speciation mechanism for the truncorotaliids was already proposed by Lazarus et al. (1995), based on biogeography, but without a definitive test.

In the modern ocean, the Globorotalia merotumida–tumida lineage is represented by G. tumida and G. ungulata, which are distributed in tropical to temperate regions. Genetic analysis has revealed that they are two distinct biological species with a single genotype each (Schiebel and Hemleben, 2017). Little is known about the ecology of these two species, although G. tumida is known to dwell in subsurface waters at the deep chlorophyll maximum (Schiebel and Hemleben, 2017). According to the phylogeny of Aze et al. (2011), Menardella menardii gave rise to the Globorotalia merotumida–tumida lineage around 9 Ma. However, Menardella menardii is absent in our late Miocene samples, while other menardellids, such as Menardella limbata and M. pertenuis, are common at several of our investigated mid-latitude to tropical sites. By 4.5 Ma, Globorotalia tumida had evolved and is common at all of our low-latitude sites (Fig. 9), while M. menardii is extremely rare, occurring only at Sites U1338 and U1490 (0.1 %), becoming more common only by the Holocene. Based on this biogeographic pattern, we propose that the G. merotumida–tumida lineage originated from a late Neogene menardellid, such as, for instance, M. limbata, which is morphologically very similar to G. merotumida. Given that G. plesiotumida and G. tumida display a deeper habitat (> 200 m) than M. limbata and other Miocene menardellids (100–200 m), we suggest that such transition may have happened through depth sympatry in tropical waters, with forms which remained reproductively isolated in the twilight zone and generated the G. merotumida–tumida lineage. Morphometric measurements on M. limbata and G. merotumida shells are required to test for an evolutionary relationship between these two species. According to our depth habitat reconstruction, G. plesitumida and G. tumida occupied a similar depth habitat at 4.5 Ma, so it is not clear from our dataset which evolutionary mechanism may have led to the origination of the latter from the former. Hull and Norris (2009) analysed speciation within this lineage and suggested that the evolution from G. plesiotumida to G. tumida happened abruptly within 44 kyr. Globorotalia ungulata appears in our record by the late Pliocene and often displays a habitat shallower than G. tumida, suggesting depth sympatry as the evolutionary mechanism leading from G. tumida to G. ungulata. However, because depth habitat reconstructions for these two species are more variable and shallower than that of other twilight zone species, more data are required to more conclusively infer speciation mechanisms.

The globoconellids originated in the late early Miocene (∼ 17 Ma) with Globoconella miozea, which is considered to descend directly from Hirsutella praescitula (Kennett and Srinivasan, 1983; Scott et al., 1990; Norris et al., 1994; Aze et al., 2011). Globoconella miotumida originated in the middle Miocene and, after the extinction of G. miozea at about 10 Ma, remained the only representative of the globoconellids until the latest Miocene. At this time, the evolutionary turnover within the group accelerated, and a number of different morphospecies originate from G. miotumida and go extinct very rapidly, until the appearance in the late Pliocene of G. inflata, which persists until today (Wei and Kennett, 1988; Wei, 1994; Aze et al., 2011). Compared to other Neogene to present taxa, the globoconellids display a more restricted geographic distribution throughout their evolutionary history. They tend to be common at mid-latitude hydrographic fronts (Schneider and Kennett, 1999; Schiebel and Hemleben, 2017; Lam and Leckie, 2020a; Brombacher et al., 2021), except for G. puncticulata and G. inflata, which extend into low-latitude regions (Norris et al., 1994). The geographic distribution of globoconellids as shown here suggests that this group was already specialized to live at hydrographic fronts in the middle Miocene, possibly feeding on phytoplankton. Starting at ∼ 5.5 Ma, cooling and the possibility to feed on sinking detritus in deeper waters (Boscolo-Galazzo et al., 2021) may have stimulated evolutionary turnover within this otherwise rather static group. The closely spaced temporal succession of morphospecies at this time may reflect ongoing evolutionary experiments in an attempt to profit from new ecological possibilities opening up at depth and outside the area of the group. The depth habitat reconstructions for G. puncticulata and G. inflata suggests that, from the Pliocene, this group started to progressively adapt to greater depths, consistent with the distinctive change in morphology between G. sphericomiozea (and other Miocene globoconellids) and its Pliocene descendants G. puncticulata and G. inflata (Kennett and Srinivasan, 1983), as well as with previous stable isotope reconstructions for these species (Schneider and Kennett, 1996). We suggest that an evolutionary transition began with the morphospecies G. puncticulata, transitional between G. sphericomiozea and G. inflata, and led to the late Pliocene speciation of G. inflata. It is not clear from our data whether depth sympatry or allopatry allowed the speciation of G. inflata, as G. puncticulata and G. inflata show similar depth habitat in our record. It may have been a combination of the two, given that G. inflata genotypes display a characteristic allopatric distribution in the ocean (Morard et al., 2011).

Our data indicate that combining stable isotopes and model-derived water column temperature is a promising approach to quantifying the depth habitat of extinct planktonic foraminiferal species. When combined with abundance and biogeographic data, depth habitat reconstructions offer insights into speciation mechanisms not resolvable with the use of one technique alone (e.g. stable isotopes). Our reconstructions indicate that both allopatry and depth sympatry played a role in the origin of modern deep-dwelling planktonic foraminiferal species. Both allopatry and depth sympatry appear to have been involved with the cladogenesis of the truncorotaliid and the globorotaliid lineages, while depth sympatry seems to be mostly involved for intra-lineage speciation.

5 Conclusions

Our global abundance and biogeographic data, combined with our depth habitat reconstructions, allow us to piece together the environmental drivers behind speciation in two of the most extensively studied groups of pelagic microfossils, planktonic foraminifera, and calcareous nannofossils over the last 15 Myr. The evolution of the new Neogene deep-water lineages of the hirsutellids, globorotaliids, globoconellids, and truncorotaliids, and of nannoplankton lower-euphotic zone and sub-euphotic zone species, resulted in the vertical stratification of species seen in the modern ocean, in particular at low latitudes. For planktonic foraminifera, such vertical stratification of species, hundreds of metres below the surface, originated through depth sympatry and the cladogenesis of new lineages via both sympatry and allopatry.

Our study places the evolution of modern plankton groups in the context of large-scale changes in ocean macroecology driven by the global climate dynamics of the past 15 Myr. The late Miocene to present evolutionary history of planktonic foraminifera and nannoplankton was linked, wherein increased efficiency of the biological pump with cooling since the middle Miocene was a shared evolutionary driver. Lower rates of organic matter remineralization in the upper part of the water column allowed the creation of new ecological niches in deep waters, by increasing food delivery and oxygen at depth for heterotrophic planktonic foraminifera and by clearing surface waters and augmenting the concentration of macronutrients at depth for nannoplankton.