Introduction
Biodiversity is multifaceted; so how should it be summarized succinctly? The ubiquitous starting point is to generate records of taxic abundance resulting in records of richness. Species are regarded by many researchers as the most intuitive unit of biology and the fundamental measure of diversity (e.g., Colwell et al. Reference Colwell, Coddington and Hawksworth1994; Purvis and Hector Reference Purvis and Hector2000; Mace et al. Reference Mace, Norris and Fitter2012; Hohenegger Reference Hohenegger2014), making species the most common currency for diversity studies. This preference for species-level studies is a result of species being argued to have independent evolutionary trajectories and histories (Purvis and Hector Reference Purvis and Hector2000) understandable to both researchers and the general public (Purvis and Hector Reference Purvis and Hector2000; Baum Reference Baum2009; Chiarucci et al. Reference Chiarucci, Bacaro and Scheiner2011; Reydon Reference Reydon, Casetta, Marques da Silva and Vecchi2019), which can aid conservation and public engagement efforts. Despite this preference, species do not represent a “silver bullet”: species show different amounts of intraspecific variation in the same ways as populations or genera. Genera also represent a biological reality (Mayr Reference Mayr1942) and share phenotypic and ecological traits (Aze et al. Reference Aze, Ezard, Purvis, Coxall, Stewart, Wade and Pearson2011). Furthermore, measuring diversity at the genus level means studies are less prone to identification error and more repeatable among different workers and the data is less prone to stochastic fluctuations that may or may not be of genuine biological interest (Hendricks et al. Reference Hendricks, Saupe, Myers, Hermsen and Allmon2014). If the goal is a unit that provides a robust summary of biodiversity change, genus-based levels provide an ecologically informative record of diversity in deep time (Hendricks et al. Reference Hendricks, Saupe, Myers, Hermsen and Allmon2014).
Following taxonomic scale, the next choice is what to measure within the sample itself. Common taxa are abundant, by definition detectable, have a broad distribution (Hannisdal et al. Reference Hannisdal, Haaga, Reitan, Diego and Liow2017), and are also likely to be influential components of ecosystems (Lennon et al. Reference Lennon, Koleff, Greenwood and Gaston2004; Gaston Reference Gaston2008; Hannisdal et al. Reference Hannisdal, Haaga, Reitan, Diego and Liow2017). Thus, it is hypothesized that common taxa contribute more to assemblage diversity than rare taxa (Lennon et al. Reference Lennon, Koleff, Greenwood and Gaston2004; Gaston Reference Gaston2008). Furthermore, for a taxon to have become common, there must have been a complex interplay of traits and environmental influences, as well as historical and spatial dynamics (Gaston Reference Gaston2008), that are replicated for establishment to be sufficiently frequent. Common taxa should therefore be able to inform our understanding of the drivers of biodiversity dynamics.
Yet, to be common is in itself rare. Very few taxa across global biodiversity are common (Gaston and Fuller Reference Gaston and Fuller2007; Gaston Reference Gaston2008; Hannisdal et al. Reference Hannisdal, Haaga, Reitan, Diego and Liow2017), so a large proportion of information is discarded when only common taxa are measured. Fluctuating abundances in rare taxa, that is, those with few individuals and geographically restricted distributions, may prove more ecologically informative, potentially acting as “ecosystem canaries” providing early warning signals for ecosystem collapse (Doncaster et al. Reference Doncaster, Chávez, Viguier, Wang, Zhang, Dong, Dearing, Langdon and Dyke2016) and insights into paleoceanographic change (Ishino and Suto Reference Ishino and Suto2020) and acting as the focus of conservation efforts (Gaston Reference Gaston2008). However, taxa vary spatially, directly influencing ecological dynamics and diversity patterns within each generation at any single location only (Patzkowsky and Holland Reference Patzkowsky and Holland2007). Consequently, what is rare in one sample or area may be common in another (Colwell Reference Colwell2009). Perhaps it is not what is rare that is important but instead what is absent from a sample, or so-called dark diversity (Pärtel et al. Reference Pärtel, Szava-Kovats and Zobel2011). On a theoretical level, far less is known about the functional role of rare taxa in their ecosystems (Lyons et al. Reference Lyons, Brigham, Traut and Schwartz2005), meaning they are easier to dismiss as unimportant and therefore to ignore (Chao et al. Reference Chao, Gotelli, Hsieh, Sander, Colwell and Ellison2014b).
To be common, rare, or absent, is a relative measure (Preston Reference Preston1948), and relative abundance requires the counting of everything to make such conclusions. Biodiversity is complex and exists on a continuum in multiple dimensions that consequently cannot be comprehensively summarized by a single number (Purvis and Hector Reference Purvis and Hector2000; Colwell Reference Colwell2009; Reich et al. Reference Reich, Tilman, Isbell, Mueller, Hobbie, Flynn and Eisenhauer2012) or straightforward categories of when a species is deemed sufficiently common to make a detectable impact on its community. Presenting diversity in integrated ways is an ideal solution (Ellison Reference Ellison2010). Such methods have been applied: effective numbers, or Hill numbers (Hill Reference Hill1973), integrate richness, evenness, and dominance in one encompassing image (Fig. 1). A drawback of effective numbers is the need for large samples of individuals. To this end, effective numbers have been applied to a range of modern-day taxonomic groups, including tropical ants (Chao et al. Reference Chao, Gotelli, Hsieh, Sander, Colwell and Ellison2014b), spiders (Chao et al. Reference Chao, Gotelli, Hsieh, Sander, Colwell and Ellison2014b, Reference Chao, Kubota, Zelený, Chiu, Li, Kusumoto, Yasuhara, Thorn, Wei, Costello and Colwell2020), corals (Chao et al. Reference Chao, Kubota, Zelený, Chiu, Li, Kusumoto, Yasuhara, Thorn, Wei, Costello and Colwell2020), and bacteria (Kang et al. Reference Kang, Rodrigues, Ng and Gentry2016). In the fossil record, assemblages are at the mercy of time and preservation (Jackson and Blois Reference Jackson and Blois2015). Therefore, abundant taxa in deep-time samples may not represent true abundance but rather a taphonomically biased sample. Despite these challenges, effective numbers have been applied meaningfully to paleoecological questions focused in the Quaternary investigating the climate and anthropogenic impacts on diversity in shallow-marine ostracods (Hong et al. Reference Hong, Yasuhara, Iwatani, Chao, Harnik and Wei2021), deep-sea ostracods (Yasuhara et al. Reference Yasuhara, Doi, Wei, Danovaro and Myhre2016), and pelagic planktonic foraminifera (Yasuhara et al. Reference Yasuhara, Wei, Kucera, Costello, Tittensor, Kiessling, Bonebrake, Tabor, Feng, Baselga, Kretschmer, Kusumoto and Kubota2020), as well as Paleozoic marine radiations (Rasmussen et al. Reference Rasmussen, Kröger, Nielsen and Colmenar2019).
Figure 1. Schematic representation of Hill numbers and how q is related to D. The color gradient in the top panel represents the weight given to abundance, with greater weight given as you move to the right. The simplified relationship between diversity line shape and the underlying assemblage is shown. Each colored dashed line is generated from an assemblage containing four taxa. The components of Assemblage 1 are represented equally within the assemblage, so the resulting diversity line is horizontal. The y-intercept is the same for both assemblages, as they have the same number of taxa (four), but Assemblage 2 has a steep gradient, as the purple taxon is more abundant than the green, red, or orange taxa. The silhouettes represent typical planktonic foraminifera of our study interval.
Conceptually, Hill numbers are the effective number of equally abundant taxa required to give the same diversity presented in the sample (Hill Reference Hill1973; Jost Reference Jost2010a; Chao et al. Reference Chao, Gotelli, Hsieh, Sander, Colwell and Ellison2014b, Reference Chao, Kubota, Zelený, Chiu, Li, Kusumoto, Yasuhara, Thorn, Wei, Costello and Colwell2020). While Hill numbers, like traditional indices such as Shannon's index (H S) and Simpson's index (H GS), can be presented as single numbers, they normally present diversity (D; Fig. 1) as a function of q, which determines how rare taxa are weighted in relation to abundant taxa (Fig. 1). Therefore, the best representation of Hill numbers is as a function of q. In uneven assemblages, this line is typically a nonlinear curve (Fig. 1) that links the three traditional indices in one image. In addition to being an integrative measure of diversity, Hill numbers also obey the replication principle (Hill Reference Hill1973). The replication principle is the requirement that when two equal assemblages with no shared taxa and equivalent relative abundances are combined, the diversity of the pooled assemblage is doubled (Hill Reference Hill1973; Chiu and Chao Reference Chiu and Chao2014). This fundamental principle is not obeyed in entropy measures such as Shannon's index. The replicable nature of Hill numbers makes them suitable for detecting diversity changes as a result of environmental perturbations whether they be anthropogenic such as oil spills (McClain et al. Reference McClain, Nunnally and Benfield2019; Miller et al. Reference Miller, Techtmann, Joyner, Mahmoudi, Fortney, Fordyce, GaraJayeva, Askerov, Cravid, Kuijper, Pelz and Hazen2020; Heritier-Robbins et al. Reference Heritier-Robbins, Karthikeyan, Hatt, Kim, Huettel, Kostka, Konstantinidis and Rodriguez-R2021) or, as in the present study, geologically transient climatic events. The commonality of units at all levels of q means that inferences can be made regarding magnitudes of change (Jost Reference Jost2007, Reference Jost2010a; Chao et al. Reference Chao, Gotelli, Hsieh, Sander, Colwell and Ellison2014b) and sample and locality differences (Hill Reference Hill1973; Chao et al. Reference Chao, Chiu and Jost2014a) and enables the transformation to commonly used general entropy metrics such as Shannon's index and Simpson's index. In addition, Hill numbers can be applied to other aspects of diversity such as phylogenetic (Chao et al. Reference Chao, Chiu and Jost2010), functional (Chiu and Chao Reference Chiu and Chao2014), and taxonomic (Chao et al. Reference Chao, Chiu and Jost2014a) diversity, with straightforward bootstrapping techniques to quantify how high proportions of singletons increase uncertainty (Chao et al. Reference Chao, Chiu and Jost2014a).
To be meaningful, Hill number calculations require sufficient and careful sampling protocols. Here we calculate Hill numbers for a deep-time community, outlining best practices for sample analysis by tracking planktonic foraminiferal diversity changes across the middle Eocene climatic optimum (MECO), ~40 Ma (Bohaty and Zachos Reference Bohaty and Zachos2003; Bohaty et al. Reference Bohaty, Zachos, Florindo and Delaney2009; Rivero-Cuesta et al. Reference Rivero-Cuesta, Westerhold, Agnini, Dallanave, Wilkens and Alegret2019; Edgar et al. Reference Edgar, Bohaty, Coxall, Bown, Batenburg, Lear and Pearson2020). The requirement for large samples of individuals means fossilized planktonic foraminifera are an ideal candidate for Hill numbers (Yasuhara et al. Reference Yasuhara, Wei, Kucera, Costello, Tittensor, Kiessling, Bonebrake, Tabor, Feng, Baselga, Kretschmer, Kusumoto and Kubota2020) as a result of their readily preserved calcium carbonate tests. In the modern oceans, planktonic foraminifera are represented by ~50 species and ~24 genera (Schiebel and Hemleben Reference Schiebel, Hemleben, Schiebel and Hemleben2001; Kucera Reference Kucera, Hillaire–Marcel and De Vernal2007; Brummer and Kučera Reference Brummer and Kučera2022), which upon death are deposited on the seafloor in vast quantities, producing ~2 Gt of calcite per year (Schiebel and Hemleben Reference Schiebel, Hemleben and Steele2008). Deposition of planktonic foraminifera has occurred nearly continuously since their evolution ~200 Ma during the Jurassic period (Fraass et al. Reference Fraass, Kelly and Peters2015), and foraminifera-rich sediments have been recovered around the globe by the coring efforts of the International Ocean Discovery Program (IODP) and its predecessors. Because planktonic foraminiferal diversity shows a strong affinity to climatic fluctuations (Ezard et al. Reference Ezard, Aze, Pearson and Purvis2011; Fraass et al. Reference Fraass, Kelly and Peters2015; Fenton et al. Reference Fenton, Pearson, Jones, Farnsworth, Lunt, Markwick and Purvis2016a; Yasuhara et al. Reference Yasuhara, Tittensor, Hillebrand and Worm2017) with a highly temporally and spatially resolved record (Fenton et al. Reference Fenton, Woodhouse, Aze, Lazarus, Renaudie, Dunhill, Young and Saupe2021), this is an ideal study system to investigate ecosystem responses to transient and rapid climatic perturbations.
Here we apply Hill numbers to understand planktonic foraminifera community response through the MECO. The middle to late Eocene encapsulates the long-term cooling from the Eocene “hothouse” of the early Eocene climatic optimum (EECO, 53–48 Ma; Westerhold et al. Reference Westerhold, Röhl, Donner and Zachos2018, Reference Westerhold, Marwan, Drury, Liebrand, Agnini, Anagnostou, Barnet, Bohaty, Vleeschouwer, Florindo, Frederichs, Hodell, Holbourn, Kroon, Lauretano, Littler, Lourens, Lyle, Pälike, Röhl, Tian, Wilkens, Wilson and Zachos2020) through to the “icehouse” of the Oligocene that started at the Eocene/Oligocene transition (EOT, 34 Ma; Westerhold et al. Reference Westerhold, Marwan, Drury, Liebrand, Agnini, Anagnostou, Barnet, Bohaty, Vleeschouwer, Florindo, Frederichs, Hodell, Holbourn, Kroon, Lauretano, Littler, Lourens, Lyle, Pälike, Röhl, Tian, Wilkens, Wilson and Zachos2020; Hutchinson et al. Reference Hutchinson, Coxall, Lunt, Steinthorsdottir, de Boer, Baatsen, von der Heydt, Huber, Kennedy-Asser, Kunzmann, Ladant, Lear, Moraweck, Pearson, Piga, Pound, Salzmann, Scher, Sijp, Śliwińska, Wilson and Zhang2021) with the establishment of continent-wide Antarctic glaciation (Zachos et al. Reference Zachos, Quinn and Salamy1996; Coxall et al. Reference Coxall, Wilson, Pälike, Lear, Backman and Arrhenius2005). The cooling trend in global temperature through this interval was interrupted by a transient (~270–500 kyr) warming event between ~40.6 and 40 Ma known as the MECO (Bohaty and Zachos Reference Bohaty and Zachos2003; Bohaty et al. Reference Bohaty, Zachos, Florindo and Delaney2009; Rivero-Cuesta et al. Reference Rivero-Cuesta, Westerhold, Agnini, Dallanave, Wilkens and Alegret2019). During the MECO, there was a transient ~3°C–6°C rise in surface and deep-water temperatures (Bohaty and Zachos Reference Bohaty and Zachos2003; Bohaty et al. Reference Bohaty, Zachos, Florindo and Delaney2009; Bijl et al. Reference Bijl, Houben, Schouten, Bohaty, Sluijs, Reichart, Damsté and Brinkhuis2010; Galazzo et al. Reference Galazzo, Thomas, Pagani, Warren, Luciani and Giusberti2014; Cramwinckel et al. Reference Cramwinckel, van der Ploeg, Bijl, Peterse, Bohaty, Röhl, Schouten, Middelburg and Sluijs2019; Henehan et al. Reference Henehan, Edgar, Foster, Penman, Hull, Greenop, Anagnostou and Pearson2020), reduced surface ocean pH (Henehan et al. Reference Henehan, Edgar, Foster, Penman, Hull, Greenop, Anagnostou and Pearson2020), and a shoaling of the calcium carbonate compensation depth (CCD; Bohaty and Zachos Reference Bohaty and Zachos2003; Bohaty et al. Reference Bohaty, Zachos, Florindo and Delaney2009). The MECO is terminated by a rapid return to pre-MECO conditions (Bohaty et al. Reference Bohaty, Zachos, Florindo and Delaney2009) and the continuation of the long-term cooling trend from the Eocene hothouse to the Oligocene icehouse.
In conjunction with the Eocene climate transition from greenhouse to icehouse conditions, there were also profound changes in planktonic foraminiferal diversity (Steineck Reference Steineck1971; Boersma and Premoli Silva Reference Boersma and Silva1986; Boersma and Silva Reference Boersma and Silva1991; Keller et al. Reference Keller, MacLeod, Barrera, Prothero and Berggren1992; Wade Reference Wade2004; Sexton et al. Reference Sexton, Wilson and Pearson2006; Wade and Pearson Reference Wade and Pearson2008; Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; Ezard et al. Reference Ezard, Aze, Pearson and Purvis2011; Galazzo et al. Reference Galazzo, Thomas, Pagani, Warren, Luciani and Giusberti2014; Fenton et al. Reference Fenton, Pearson, Jones, Farnsworth, Lunt, Markwick and Purvis2016a). Middle Eocene biotic changes in planktonic foraminifera include: (1) the progressive extinction of surface-dwelling, symbiont-bearing taxa (Boersma and Premoli Silva Reference Boersma and Silva1986; Boersma and Silva Reference Boersma and Silva1991; Keller et al. Reference Keller, MacLeod, Barrera, Prothero and Berggren1992; Wade Reference Wade2004; Wade and Pearson Reference Wade and Pearson2008); (2) a reduction in test size (Schmidt et al. Reference Schmidt, Thierstein and Bollmann2004; Wade and Pearson Reference Wade and Pearson2008; Wade and Olsson Reference Wade and Olsson2009); (3) development of latitudinal size (Schmidt et al. Reference Schmidt, Thierstein and Bollmann2004) and diversity (Fenton et al. Reference Fenton, Pearson, Jones, Farnsworth, Lunt, Markwick and Purvis2016a) gradients alongside major assemblage fluctuations (Steineck Reference Steineck1971; Keller Reference Keller1983; Boersma and Premoli Silva Reference Boersma and Silva1986; Boersma et al. Reference Boersma, Silva and Shackleton1987; Hallock et al. Reference Hallock, Silva and Boersma1991; Keller et al. Reference Keller, MacLeod, Barrera, Prothero and Berggren1992; Sexton et al. Reference Sexton, Wilson and Pearson2006; Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; Galazzo et al. Reference Galazzo, Thomas, Pagani, Warren, Luciani and Giusberti2014); and (4) changes in ecology, for example, loss or inhibition of algal photosymbionts from hosting taxa (Wade et al. Reference Wade, Al-Sabouni, Hemleben and Kroon2008; Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013) and shallowing depth habitat of Hantkenina (Coxall et al. Reference Coxall, Pearson, Shackleton and Hall2000). Yet our understanding of planktonic foraminifera ecosystem dynamics across the MECO remain relatively understudied compared with other periods of the Eocene such as the EOT (e.g., Pearson et al. Reference Pearson, Mcmillan, Wade, Jones, Coxall, Bown and Lear2008; Wade and Pearson Reference Wade and Pearson2008; Pearson and Wade Reference Pearson and Wade2015). The MECO resulted in a global crisis for muricate taxa (Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013), varying symbiotic taxon responses (Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013; Gebhardt et al. Reference Gebhardt, Ćorić, Darga, Briguglio, Schenk, Werner, Andersen and Sames2013; Arimoto et al. Reference Arimoto, Nishi, Kuroyanagi, Takashima, Matsui and Ikehara2020; Kearns et al. Reference Kearns, Bohaty, Edgar, Nogué and Ezard2021), and increased abundance of ecologically flexible (Galazzo et al. Reference Galazzo, Thomas, Luciani, Giusberti, Frontalini and Coccioni2015; Kearns et al. Reference Kearns, Bohaty, Edgar, Nogué and Ezard2021) and small opportunistic taxa (Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010). What we lack, however, is an integrated assemblage perspective on these idiosyncratic changes, pieced together from different sampling localities. Here, using Hill numbers, we generate the first midlatitude diversity record of planktonic foraminifera at North Atlantic sites through the MECO to investigate how planktonic foraminifera communities responded to the MECO and how this event may have influenced subsequent extinction events observed in the late Eocene. Furthermore, we analyze diversity within two size fractions (>63 μm and >180 μm) to understand the effects of sampling bias on diversity and its implications for our understanding of biotic responses to climatic perturbations.
Materials and Methods
Materials
IODP Expedition 342 targeted clay-rich Paleogene sediment drifts ~700 km east-southeast of Newfoundland in the northwest Atlantic Ocean (Norris et al. Reference Norris, Wilson and Blum2014), which were deposited at a paleolatitude of ~32.5°N (Supplementary Fig. 1). Expedition 342 Sites U1406 (40°21.0′N, 51°39.0′W; modern water depth: ~3814 m), U1408 (41°26.3′N, 49°47.1′W; modern water depth: ~3022 m), and U1410 (41°19.6993′N, 49°10.1847′W; modern water depth: ~3387 m) recovered clay-rich nannofossil ooze drift deposits well above the late Paleogene CCD, providing near-continuous records of well-preserved microfossils from ~47 Ma through the Eocene and into the Oligocene (Norris et al. Reference Norris, Wilson and Blum2014; Boyle et al. Reference Boyle, Romans, Tucholke, Norris, Swift and Sexton2017). Using low-resolution bulk stable isotope data (S. M. Bohaty unpublished data) as a guide, cores from Sites U1406, U1408, and U1410 were sampled to capture a 7 Myr interval of the middle Eocene spanning the MECO. In total, 33 samples of 25 cc between 38 and 45 Ma were studied. Due to changes in sediment accumulation rates during parts of the MECO, sampling resolution ranges from ~20 kyr during the MECO to ~900 kyr outside the MECO.
Sample ages from Sites U1408 and U1410 were calculated based on age–depth models constructed using available biostratigraphy and magnetostratigraphy (Norris et al. Reference Norris, Wilson and Blum2014). The 2012 geological timescale was then used for age calibrations for the middle Eocene geomagnetic reversals (GTS2012; Gradstein et al. Reference Gradstein, Ogg, Schmitz and Ogg2012). Sample ages for Site U1406 are based upon shipboard biostratigraphic and magnetostratigraphic data (Norris et al. Reference Norris, Wilson and Blum2014; Van Peer Reference Van Peer, Liebrand, Xuan, Lippert, Agnini, Blum, Blum, Bohaty, Bown, Greenop and Kordesch2017). Sample information, including calculated ages, is presented in Supplementary Table 1.
Sample Preparation
The sample material was disaggregated in a sodium hexametaphosphate solution and then washed over a 38 μm sieve with milli-Q water until the water ran clear. Following 24 hours of drying in a low-temperature oven (<50°C), samples were weighed to determine the weight percent coarse fraction (>38 μm). Subsequently, each sample was split, using a microsplitter, providing two representative halves: one for diversity analysis (this study) and the other for geochemical analysis (Kearns et al. Reference Kearns, Bohaty, Edgar, Nogué and Ezard2021). The half reserved for diversity analysis in this study was then split again to allow analyses at two different size fractions. Planktonic foraminiferal assemblage studies typically only analyze size fractions >150 μm (Kucera et al. Reference Kucera, Weinelt, Kiefer, Pflaumann, Hayes, Weinelt, Chen, Mix, Barrows, Cortijo, Duprat, Juggins and Waelbroeck2005) to avoid sampling juvenile specimens and to enable species-level identification (Al-Sabouni et al. Reference Al-Sabouni, Kucera and Schmidt2007, Reference Al-Sabouni, Fenton, Telford and Kučera2018). This, by definition, biases assemblages toward larger forms, despite suggestions that analyzing a >63 μm size fraction, especially in polar regions where species are generally smaller, is more representative of true diversity (Al-Sabouni et al. Reference Al-Sabouni, Kucera and Schmidt2007). To test whether a smaller size fraction is more characteristic of diversity at midlatitude, nonpolar sites like IODP Expedition 342, we determined diversity in two size fractions: >63 μm and >180 μm. To avoid juveniles in the smaller size fractions, only individuals showing adult characteristics related to aperture position, keels, and fully developed pore structure in macroperforate forms were picked for analysis (Brummer et al. Reference Brummer, Hemlebent and Spindlert1986).
Diversity Analysis
A sample of 300 individuals is considered sufficient to estimate diversity in foraminiferal assemblages (Al-Sabouni et al. Reference Al-Sabouni, Kucera and Schmidt2007), despite the potential for missing rare specimens due to low abundances (Jost Reference Jost2010b). For this study, each sample in both size fractions (>63 μm and >180 μm) was further split using a microsplitter until approximately 300 individuals were present on the picking tray, with a minimum cutoff of 200 specimens. All individuals in the subsample were then picked to avoid bias as a result of uneven distribution on the tray and identified to the genus level (Supplementary Tables 2–4) based on published taxonomy (Pearson et al. Reference Pearson, Olsson, Huber, Hemleben and Berggren2006; Wade et al. Reference Wade, Olsson, Pearson, Huber and Berggren2018).
To understand diversity changes further, we then classified all genera into morphogroups (Supplementary Tables 5, 6) adapted from previous classifications (Aze et al. Reference Aze, Ezard, Purvis, Coxall, Stewart, Wade and Pearson2011) and depth habitats (Supplementary Tables 7, 8). We based morphogroup classifications on morphological traits (Supplementary Table 9) and depth habitats (Supplementary Table 10) on published ecological inferences obtained from stable isotope measurements (summarized in Pearson et al. Reference Pearson, Olsson, Huber, Hemleben and Berggren2006; Wade et al. Reference Wade, Olsson, Pearson, Huber and Berggren2018). Relative abundances and effective diversity curves were than calculated for each genus, morphogroup, and ecogroup.
We calculate diversity as a curve using Hill numbers (Hill Reference Hill1973):
$${}_{}^q D = \left({\mathop \sum \limits_{i = 1}^S p_i^q } \right)^{1/( {1-q} ) }$$where S is the number of taxa and p i the frequency of the ith taxa. The value of D is dependent on the order, q, which determines how rarity is weighted in relation to abundance. At q = 0, taxic richness is measured such that abundance is ignored, as rare taxa are weighted more heavily than common taxa compared with higher powers of q (Fig. 1). As q gets larger, the weighting toward rare taxa is reduced and relative abundance is considered. At q = 1, rare and common taxa are equally weighted, which equates to the exponential of Shannon's index (Chao and Jost Reference Chao and Jost2012; Fig. 1). At q = 2, only relative abundance is accounted for, removing the influence of rare taxa, so this measure is equivalent to the inverse of Simpson's index (Chao and Jost Reference Chao and Jost2012; Fig. 1). While these integer values are useful reference points, the strength of the Hill number approach is how the continuum of q values (the slope of the effective diversity curve) can be used to understand the evenness of the assemblage. If an assemblage is made up of equal numbers of represented taxa, then the diversity curve will be flat, as abundance does not vary among between groups and no taxon is rare (Fig. 1). In contrast, if the curve has a high gradient and plummets into a plateau, then the assemblage can be interpreted as uneven with lots of rare taxa and a few dominant groups (Fig. 1).
We outline the workflow for calculating our diversity curves, which follows Chao and Jost (Reference Chao and Jost2015), in the Supplementary Material. Effective diversity was calculated at the default 0.1 intervals for q between 0 and 2 (Chao and Jost Reference Chao and Jost2015; Supplementary Table 11). Confidence intervals were generated at the 95% level for each diversity curve by bootstrapping 1000 times.
Fragmentation
A challenge to using paleoecological data is the inevitable influence of taphonomic bias. Assemblage data of planktonic foraminifera can be heavily influenced by the physiochemical process of dissolution as a result of their shell (test) composition (Berger Reference Berger1971; Malmgren Reference Malmgren1987; Nguyen et al. Reference Nguyen, Petrizzo and Speijer2009). The susceptibility of foraminifera to dissolution is strongly species specific based on the physical structure of the test wall (e.g., relative porosity and thickness; Nguyen et al. Reference Nguyen, Petrizzo and Speijer2009, Reference Nguyen, Petrizzo, Stassen and Speijer2011; Nguyen and Speijer Reference Nguyen and Speijer2014) and is influenced by the microenvironment of the individual, which influences test chemistry and causes interspecific differences in dissolution susceptibility (Berger Reference Berger1970; Nguyen et al. Reference Nguyen, Petrizzo, Stassen and Speijer2011; Petro et al. Reference Petro, Pivel and Coimbra2018). To account for this variability, we use an accepted fragmentation proxy to estimate the dissolution levels (Le and Shackleton Reference Le and Shackleton1992), using the proportion of planktonic foraminiferal test fragments (Frag) and whole specimens:
$$ \eqalign{\!\!{\rm Fragmentation\;}( {\rm \% } ) = & \,\, ( {{\rm Frag}/8} ) / \cr& \,\, ( {{\rm Frag}/8\;\times {\rm Whole}} ) \times 100$$We classify a fragment as anything <75% of a whole specimen (which is more conservative than the <50% previously used; Malmgren Reference Malmgren1987). Foraminifera have a tendency to break into multiple pieces; therefore, the percentage of fragments in a sample varies nonlinearly with dissolution (Le and Shackleton Reference Le and Shackleton1992). To account for this, a divisor is used to represent the average number of pieces a foraminifera breaks into, and we follow previous work and set the divisor as 8 (Le and Shackleton Reference Le and Shackleton1992; Leon-Rodriguez and Dickens Reference Leon-Rodriguez and Dickens2010). We use a baseline of 20% fragmentation to indicate normal levels of fragmentation and dissolution (Pfuhl and Shackleton Reference Pfuhl and Shackleton2004). Samples sieved at 63 μm are expected to have higher fragmentation than samples sieved at a larger size fraction, as fragments progressively break into smaller pieces and smaller individuals are less robust. Therefore, we use the percentage of coarse fraction after sieving to assess potential dissolution effects on the assemblage, as dissolution reduces the absolute abundance of planktonic foraminifera in a sample while ecological change causes taxon relative abundance fluctuations. Fragmentation was calculated twice on 18 samples (9 samples from the >180 μm fraction and 9 samples from the >63 μm fraction) with high repeatability (92%; Supplementary Table 12).
Statistical Methods
Generalized Additive Models
Diversity has a nonlinear relationship with time. To assess the impact of sample age and size fraction on diversity, we applied nonparametric generalized additive models (GAMs) using the R package mgcv (v. 1.8.33; Wood Reference Wood2017) in the R environment (v. 4.0.3; R Core Team 2020). Before model fitting, integer values of Hill numbers were back transformed to genus richness, Shannon's index (H S), and Simpson's index (H GS) and used as response variables. Models were constructed with a smooth (nonparametric, nonlinear) function of age and a linear predictor of fragmentation to control for the impact on dissolution on diversity. Models were fit using a Gaussian distribution with an identity link using a generalized cross-validation model (GCV) method. A GCV method was used instead of restricted maximum-likelihood estimation (REML) due to the small number of samples through time (Wood Reference Wood2011). The code to obtain predictions based on observed and defined fragmentation can be found in the Supplementary Material.
Model selection among a relevant model set including a null model (Table 1) was based on Akaike information criterion corrected for small sample size (AICc) and diagnostic plots. The Supplementary Material provides further detail on back transformation, models fit (including all annotated code), and model selection information. Model results relating to significance of smoothing parameters are presented with effective degrees of freedom (edf), F-statistics (F), and p-values (p). The edf indicates the complexity of the curve; an edf of 1 indicates a straight line, while an edf of 2 indicates a quadratic curve, and so on. In addition, where appropriate, parametric coefficients are presented with the coefficient (β), t-value (t), standard error (SE), and p-value (p).
Table 1. Table showing the structure of all models fit. Diversity is replaced by genus and morphogroup for each set of models, and Akaike information criterion (AIC) weights are presented. The smooth term is denoted by s().*Best-fitting model. †Null model. All other statistical output, including df and AIC, is provided in Supplementary Tables 13–18.
Kruskal-Wallis and Dunn Tests
To investigate differences in Hill numbers in response to paleoclimatic and paleoceanographic changes, samples were divided into three time slices representing different climate phases (pre-MECO: >41.94 Ma; MECO: 41.09–40.14 Ma; post-MECO: <40.14 Ma). The difference between these intervals was assessed when q <1 (weighted toward rarity) and q >1 (relative abundance taken into account) using a Kruskal-Wallis test. A Kruskal-Wallis test was used to investigate whether a difference was present between intervals, and additionally the effect size of the intervals was calculated based on the H statistic from the Kruskal-Wallis test. Following detection of a statistically significant impact of interval in the Kruskal-Wallis test (p < 0.01), a post hoc Dunn test using the R package “FIA” (Ogle et al. Reference Ogle, Doll, Wheeler and Dinno2022) was applied, due to unequal observations in each interval (Zar Reference Zar2010), to identify intervals which were significantly different from each other.
Results
Fragmentation
The degree of fragmentation varies across our record from 1.34% to 30.80% (Supplementary Fig. 2, Supplementary Table 11), with generally increased fragmentation in the smaller size fraction, as expected. In total, seven samples were above the baseline “normal” fragmentation of 20% (Pfuhl and Shackleton Reference Pfuhl and Shackleton2004), of which six were in the >63 μm size fraction. These were primarily within the MECO interval (~41‒40 Ma) and at ~38 Ma.
Traditional Diversity Indices
In total 18 genera consisting of 11 morphogroups occupying three depth habitats were identified (Supplementary Table 2). On average, samples consisted of 11 and 9 genera in the >63 μm and >180 size fractions, respectively (Supplementary Tables 3, 4). In the >63 μm fraction, only three genera were present in all 33 samples (Subbotina, Acarinina, and Planorotalites; Supplementary Table 3), while in the >180 μm size fraction, only Subbotina was found in all samples (Supplementary Table 4). This study is based on genera, for reasons outlined in the “Introduction”; however, we note that all genera were represented by approximately two or fewer species, except for Subbotina, which was represented by approximately six species.
Integers of effective diversity are equivalent to transformed versions of common diversity measures (genus richness [q = 0], Shannon's index (H S) [q = 1], and Simpson's index (H GS) [q = 2]). To understand how commonly used diversity indices changed through time and aid comparison with other studies, we back-transformed calculated Hill numbers into genus richness, H S and H GS indices for genera, morphogroup, and ecogroup, and fitted GAMs (Fig. 2, Table 1). For all diversity indices, based on AICc, the best-fitting GAMs (Tables 1, 2, Supplementary Tables 13–21, Supplementary Figs. 3–11) suggest a change in diversity as a function of size fraction, with varying intercepts for size fraction (with the smaller size fraction giving consistently higher values), during the Eocene for richness, H S, and H GS. We concentrate here on generic and morphogroup changes in terms of richness, H S, and H GS (depth-habitat effects on diversity changes are discussed in “Hill Numbers and Relative Abundance Fluctuations”). Depth-habitat analysis showed similar patterns, but with only three depth groups, the changes observed were nonconsequential in terms of H S and H GS and are therefore provided in Supplementary Figs. 9–12 and Supplementary Tables 21–23. The ΔAICc between the best-fitting models for genera and morphogroup and the next best models ranged from 11.595 to 71.331 (Supplementary Tables 14–18), implying that the second-ranked models had “essentially no” support (Burnham and Anderson Reference Burnham and Anderson2002).
Table 2. Table showing the parametric coefficients for best-fitting generalized additive models (GAMs) shown in Table 1 based on Akaike information criterion (AIC) difference compared with next best fitting model being more than 2. Frag, fragmentation; size, size fraction (>63 μm or >180 μm).
Figure 2. North Atlantic predicted diversity from International Ocean Discovery Program (IODP) Expedition 342 (Sites U1406, U1408, and U1410) as a function of time at >63 μm (green) and >180 μm (purple) size fractions. Note that these predictions control for fragmentation, which we assume in this prediction to be 10% (the mean fragmentation of our samples). Raw data are shown as filled circles. Black shaded area represents 95% confidence intervals around the central predicted response. A–C show generic diversity indices; D–F show morphological diversity indices. The light gray box represents the middle Eocene climatic optimum (MECO) interval.
Following model selection, we used the best-fitting model (Tables 1, 2) to predict diversity values across our study interval at a mean fragmentation of 10% based on the mean of our fragmentation counts (10.37%; Supplementary Fig. 2, Supplementary Table 11) to produce diversity curves and 95% confidence intervals (Fig. 2). For completeness, model predictions were also done at 5% and 20% fragmentation, which resulted in no change in predicted diversity values.
Richness
Spline complexity (“wigglyness”) for genera differed between size fractions, with a more complex spline predicted for the >63 μm size fraction (edf = 8.54, F = 3.10, p < 0.01) compared with the >180 μm size fraction (edf = 2.65, F = 4.27, p < 0.01) (Fig. 2A). A similar pattern was observed in the morphogroup models, where the predicted spline for >63 μm size fraction is more complex (edf = 8.26, F = 2.16, p < 0.05) than the >180 μm size fraction (edf = 2.38, F = 3.10, p < 0.05) (Fig. 2D). The complex nature of the >63 μm size fraction spline illustrates intersample variability represented in the larger confidence intervals compared with the >180 μm size fraction (Fig. 2A,D).
Morphological and generic richness profiles generally follow a similar pattern, with increasing richness initially between 45 and 44 Ma, followed by a period of relative stasis until ~41.5 Ma (Fig. 2A,D). In the >63 μm size fraction, generic and morphological richness peaked at ~40.55 Ma, coinciding with the early stages of the MECO (generic: 14.71 ± 0.98, morphological: 8.45 ± 0.55; Fig 2A,D). In the >180 μm size fraction, the peak in richness is much less pronounced and occurs ~1 Myr before the MECO at 41.54 Ma (generic: 9.52 ± 0.47; Fig. 2A) and 41.89 Ma (morphological: 6.14 ± 0.30; Fig. 2D). Peaks in morphological and generic richness in the >180 μm size fraction are followed by a decline of a mean of 1.64 morphogroups and 2.06 genera by the end of our record at 38.00 Ma (Fig. 2A,D). The wide 95% confidence intervals around the richness declines (Fig. 2A,D) suggest no detectable fall, as the confidence intervals could encapsulate a straight horizontal line. In contrast, the >63 μm size fraction shows a greater degree of intrasample variability resulting in more complex GAMs that predict a large decline in morphological (−2.17 morphogroups) and generic (−5.62 genera) richness following the MECO at ~40.5 Ma (Fig. 2A,D).
The most influential predictor of generic and morphological richness was size, with a predicted reduction in overall richness of 2.537 genera (β = −2.537, SE = 0.468, t = −5.423, p < 0.001) and 1.728 morphogroups (β = −1.728 SE = 0.277, t = −6.235, p < 0.001), calculated through analysis assemblages from the >180 μm size fraction rather than >63 μm size fraction (Table 2). This means that 2.537 genera and 1.728 morphogroups represented in the >63 μm size fraction are not present in the >180 μm size fraction. Fragmentation was also a significant predictor for generic richness with a predicted 0.14 decrease in richness per 1% increase in fragmentation (β = −0.14, SE = 0.048, t = −13.16, p < 0.001; Table 2).
Shannon's Index
The predicted curves for H S are smoother than those for richness (Fig. 2). However, unlike richness, the model predicted a more complex age spline for generic H S in the >180 μm size fraction (edf = 3.06, F = 11.65, p < 0.001) than in the >63 μm size fraction (edf = 2.18, F = 2.73, p > 0.05; Fig. 2B, Supplementary Table 20). Among genera, size is the only significant predictor of diversity (β = −0.58, SE = 0.10, t = −8.68, p < 0.001; Supplementary Table 20), with the >180 μm size fraction predicted to increase to a peak of 1.64 ± 0.08 at 41.89 Ma followed by a steep decline until 38.00 Ma (Fig. 2B). In contrast, the >63 μm size fraction gradually increases, reaching a maximum H S of 1.99 ± 0.07 at 42.10 Ma (Fig. 2B).
For morphological H S, the age spline for the >63 μm fraction does not differ detectably from a straight line (edf = 1.62, F = 2.28, p > 0.05; Supplementary Table 20), in contrast to the wiggly spline for >180 μm (edf = 3.47, F = 13.67, p < 0.001; Fig. 2E, Supplementary Table 20). Both fragmentation and size fraction are significant predictors (fragmentation: β = −0.011, SE = 0.004, t = −2.74, p < 0.01; size fraction: β = −0.61, SE = 0.05, t = −11.61, p < 0.001; Table 2), but size fraction has a larger, more meaningful impact on diversity, with a reduction of 0.61 morphogroups in the >180 μm size fraction compared with the >63 μm size fraction. The peak in >180 μm morphological H S (1.18 ± 0.07) is predicted at 42.40 Ma, 0.50 Ma before the peak in generic H S in the same size fraction (Fig. 2B,E).
Simpson's Index
Through the narrow range of values allowed for available values for H GS (between 0 and 1), intersample variation was high (Fig. 2 C,F), and the predicted spline follows a pattern similar to that in H S (Fig. 2). The GAMs predicted a complex age spline for genera (edf = 3.31, F = 19.11, p < 0.001; Supplementary Table 20) and morphogroup (edf = 3.64, F = 17.94, p < 0.001; Supplementary Table 20) H GS in the >180 μm size fraction. Both splines reach peaks before the MECO at 41.82 Ma (Fig. 2C) and 42.31 Ma (Fig. 2F) for morphogroup and generic H GS, respectively. The age spline for the smaller size fractions is not detectably different from straight lines (genera: edf = 1.88, F = 2.01, p > 0.05; morphogroup: edf = 1.06, F = 2.24, p > 0.05; Fig. 2C,E, Supplementary Table 20). In genera, only fragmentation is a significant predictor equating to a predicted 0.21 ± 0.02 reduction in generic Simpson's index per 1% increase in fragmentation (β = −0.21, SE = 0.02, t = −9.40, p < 0.001; Table 2), while for morphogroup, both fragmentation (β = −0.01, SE = 0.002, t = −3.06, p < 0.01; Table 2) and size fraction (β = −0.28, SE = 0.02, t = −11.62, p < 0.001; Table 2) are significant predictors of H GS.
Hill Numbers and Relative Abundance Fluctuations
Genera
When relative abundance is not considered (q < 1), pre-MECO (>41.09 Ma), MECO (41.09–40.14 Ma), and post-MECO (<40.14 Ma) intervals are all different from each other (Fig. 3A,D; p < 0.001), but the effect of size (η2[H] = 0.18) and magnitude of differences between intervals (pre-MECO, MECO, and post-MECO) is only large (defined in the statistical test as an effect size (η2[H]) > 0.14) in the >63 μm size fraction (Supplementary Table 24). This suggests substantial differences in absolute numbers of genera through the middle Eocene, with highest values in the MECO (41.09–40.14 Ma) followed by a decline into post-MECO (<40.14 Ma) assemblages below pre-event levels (Fig. 3A,D).
Figure 3. North Atlantic diversity curves showing Hill number calculations based on abundance counts presented in Supplementary Tables 3 and 4. Morphogroup and depth habitat follow the classification outlined in Supplementary Tables 9 and 10. A–C reflect diversity changes at the >63 μm size fraction, while D–F reflect changes at the >180 μm size fraction. All panels show a reduction in evenness in post–middle Eocene climatic optimum (post-MECO) communities compared with pre-MECO and MECO assemblages. D–F, The larger size fraction has steeper curves, reiterating the potential dangers of unrepresentative community sampling. Vertical dotted lines are present where q = 1 and q = 2, as these correlate to the exponential of Shannon's index (q = 1) and the inverse of Simpson's index (q = 2) presented in Fig. 2. Lines are colored to represent paleoceanographic intervals. Note one horizontal yellow line in F illustrating a perfectly even assemblage. The gray bands represent 95% confidence intervals.
When relative generic abundance is considered (q >1), the effect size of time interval is only large (η2[H] = 0.25) and significant (p < 0.01) in the >180 μm size fraction (Supplementary Table 23). A Dunn's test shows that a significant difference exists between the post-MECO interval (<40.14 Ma) and the other intervals (p < 0.01; Supplementary Table 24), with the MECO (41.09–40.14 Ma) and pre-MECO interval (>41.09 Ma) showing no significant differences.
Morphogroup
In the >63 μm size fraction, the effective diversity curves show only subtle separation between assemblages because of paleoceanographic changes (interval) (Fig. 3B). A Kruskal-Wallis test revealed morphological effective diversity was only significantly different (p < 0.01) between time intervals, with a large effect size of interval (η2[H] = 0.25) when relative abundance was considered (q >1; Fig. 3B, Supplementary Table 24). Based on a Dunn's test, the detectable interval differences are between the MECO (41.09–40.14 Ma) and pre-MECO (>41.09 Ma; p < 0.05), as well as between the post-MECO (<40.14 Ma) and pre-MECO (>41.09 Ma; p < 0.05; Supplementary Table 25).
In the >180 μm size fraction, the effect size of time interval is clear when rare morphologies are influential (q < 1) and when they are discounted (q = 1–2) (Supplementary Table 24). A Dunn test showed that there is no detectable difference between pre-MECO and MECO (41.09–40.14 Ma) samples (p > 0.1) at any level of q (Supplementary Table 25), resulting in overlapping effective diversity curves, except for those assemblages grouped in the post-MECO (<40.14 Ma) assemblages colored purple (Fig. 3E). Additionally, the post-MECO (<40.14 Ma) assemblages show a decrease in evenness compared with the preceding intervals (Fig. 3E).
Depth Habitat
Compared with the generic and morphogroup analyses, effective depth-habitat richness (q = 0) in assemblages is the same in both size fractions (3; Fig. 3C,F), with no differences in depth habitat as a function of size fraction or sample age. A Kruskal-Wallis test showed no difference in effective diversity between paleoceanographic intervals in the >63 μm size fraction, which is illustrated by the overlap of effective diversity curves (Fig. 3C). This implies that there was no change in depth-habitat evenness through the middle Eocene, meaning no organisms of a certain depth habitat were dominating assemblages.
In comparison, the depth-habitat effective number curves show separation in the post-MECO (<40.14 Ma) samples of the >180 μm size fraction (Fig. 3F). A Kruskal-Wallis test showed that the interval had a large effect size (η2[H] = 0.262) with a clear impact (p < 0.01) when q is between 1 and 2 in the >180 μm size fraction (Supplementary Table 24). A Dunn test show that the detectable differences are between the post-MECO (<40.14 Ma) interval and both preceding intervals (p < 0.01; Supplementary Table 25). The gradient change of the effective diversity curves also shows that the post-MECO (<40.14 Ma) samples are uneven compared with the other intervals.
Discussion
Understanding biodiversity responses to climate change is challenging, particularly in deep time. A focus on relative abundance changes and biogeographic comparisons can complicate broader interpretations because of the idiosyncratic responses of taxa to environmental change. By using Hill numbers, we have been able to generalize and assess the biodiversity response of planktonic foraminifera temporally to transient warming (Fig. 3) at midlatitudes, while maintaining the ability to investigate more specific biodiversity measures such as richness (Fig. 2) and relative abundance changes (Fig. 4).
Figure 4. Relative abundance plots of genera across the North Atlantic middle Eocene from International Ocean Discovery Program (IODP) Expedition 342 Sites U1406, U1408, and U1410 separated by depth habitat and size fraction. A, Surface dwellers with symbiont-bearing taxa filled with a pattern; B, thermocline dwellers; C, subthermocline dwellers. Note different color schemes are used per depth habitat for ease of viewing. Symbiont-bearing foraminifera: Acarinina, Morozovelloides, and Globigerinatheka are indicated by a striped pattern. The gray horizontal box represents the middle Eocene climatic optimum (MECO) interval. Turborotalita, Orbulinoides, and Hantkenina are not included in this plot, as they occurred in such low numbers (1–5 absolute abundance).
Using this approach, we show increases in morphological and generic richness coincident with the early and middle stages of MECO warming in the >63 μm size fraction (Fig. 2A,D), which is reflected in our count data with ~70% of all genera and morphogroups we found present at ~40.5 Ma (Supplementary Tables 2–8). In addition, we found that analytical choice of size fraction resulted in apparent differences in planktonic foraminifera response to both MECO warming and post-MECO cooling (Figs. 2, 3). The loss of symbiont-bearing foraminifera only changes depth-habitat diversity in the >180 μm size fraction (Fig. 3F), because these genera are replaced in the mixed layer by increased numbers of nonsymbiotic Chiloguembelina and Planorotalites in the >63 μm size fraction (Fig. 4).
Influence of Dissolution on Diversity Analysis
Dissolution has the potential to shift planktonic foraminiferal assemblages from representing environmentally shaped life assemblages to taphonomically shaped death assemblages (Berger Reference Berger1971; Thunell Reference Thunell1976), biasing climatic and biotic interpretations of these assemblages (Berger Reference Berger1973). Dissolution can be morphologically selective (Berger Reference Berger1970; Boltovskoy and Totah Reference Boltovskoy and Totah1992; Petrizzo et al. Reference Petrizzo, Leoni, Speijer, De Bernardi and Felletti2008; Nguyen et al. Reference Nguyen, Petrizzo and Speijer2009, Reference Nguyen, Petrizzo, Stassen and Speijer2011, but see Petro et al. Reference Petro, Pivel and Coimbra2018) with species-specific tendencies (Nguyen et al. Reference Nguyen, Petrizzo, Stassen and Speijer2011; but see Berger Reference Berger1970; Malmgren Reference Malmgren1987). Across the MECO, extensive dissolution as a result of shoaling CCD has been recorded in the Pacific, Indian, and Atlantic Oceans (Lyle et al. Reference Lyle, Olivarez Lyle, Backman and Tripati2005; Bohaty et al. Reference Bohaty, Zachos, Florindo and Delaney2009; Pälike et al. Reference Pälike, Lyle, Nishi, Raffi, Ridgwell, Gamage, Klaus, Acton, Anderson, Backman, Baldauf, Beltran, Bohaty, Bown, Busch, Channell, Chun, Delaney, Dewangan, Dunkley Jones, Edgar, Evans, Fitch, Foster, Gussone, Hasegawa, Hathorne, Hayashi, Herrle, Holbourn, Hovan, Hyeong, Iijima, Ito, Kamikuri, Kimoto, Kuroda, Leon-Rodriguez, Malinverno, Moore, Murphy, Murphy, Nakamura, Ogane, Ohneiser, Richter, Robinson, Rohling, Romero, Sawada, Scher, Schneider, Sluijs, Takata, Tian, Tsujimoto, Wade, Westerhold, Wilkens, Williams, Wilson, Yamamoto, Yamamoto, Yamazaki and Zeebe2012). Despite our study sites sitting well above the known late Paleogene CCD (Norris et al. Reference Norris, Wilson and Blum2014), we still find that ~11% of our samples (7 out of 66) had higher fragmentation than what is considered normal for a well-preserved sample (Pfuhl and Shackleton Reference Pfuhl and Shackleton2004). While this indicates some degree of dissolution, probably reflecting increased carbonate dissolution due to the shoaling of the lysocline during the MECO (Boscolo Galazzo et al. Reference Boscolo Galazzo, Giusberti, Luciani and Thomas2013; Savian et al. Reference Savian, Jovane, Frontalini, Trindade, Coccioni, Bohaty, Wilson, Florindo, Roberts, Catanzariti and Iacoviello2014), there was no observable drop in overall planktonic foraminifera abundances, which would indicate dissolution impacted assemblages (Malmgren Reference Malmgren1987). We do detect a statistically significant effect of dissolution in our statistical models, and recommend accounting for it in statistical analyses, but the small effect sizes (e.g., up to 55 times smaller than the effect of size fraction choice) and unchanged predictions when fragmentation is doubled and halved (Supplementary Figs. 13, 14) suggest that dissolution is not a strong driver of the diversity dynamics we report (Figs. 2–4).
Transient Climate Impacts on Specialist Feeding Ecologies
Our samples are the first midlatitude open-ocean samples analyzed for assemblages across the middle Eocene. Therefore, our results give a unique insight into the impacts of the MECO on symbiont-bearing foraminifera (Acarinina, Morozovelloides, and Globigerinatheka) and motivation for further studies at high-latitude sites.
We observe abundance decreases on the Newfoundland margin before the MECO at ~40.50 Ma in the >180 μm size fraction and at ~41.31 Ma in the >63 μm size fraction that persist post-MECO (Fig. 4, Supplementary Tables 3, 4). At the lower latitude Ocean Drilling Program Site 1051 in the North Atlantic Ocean (~25°N, Blake Nose), large Acarinina (>300 μm) abundance only temporarily decreases during peak warming of the MECO (Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013); in the subtropical Alano section, the abundance of Acarinina is high before and during the MECO, but then abruptly decreases post-MECO and remains low (Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010). Our observed decline in abundance is notably smaller on the Newfoundland margin in the North Atlantic (~20% reduction in Acarinina) compared with the subtropical Alano section. Acarinina relative abundance never recovers following the decline in our record, instead staying consistently low as in the Tethys (Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; Fig. 4), unlike the lower-latitude Blake Nose (~25°N), where Acarinina recovers in both abundance and test size (Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013). Though not as abundant as Acarinina, Morozovelloides is present in our samples, with a peak relative abundance in both size fractions at 41.31 Ma before the MECO (~20% in the >180 μm fraction and ~11% in the >63 μm fraction; Fig. 4), followed by a decline in relative abundance through MECO warming, despite being thermophilic, and to the end of our record (Fig. 4). The general trend of post-MECO reduction in relative abundance of Morozovelloides is observed at other localities (Wade et al. Reference Wade, Al-Sabouni, Hemleben and Kroon2008; Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013), although in the Tethys, Morozovelloides is scarcely abundant throughout the middle Eocene (Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; Gebhardt et al. Reference Gebhardt, Ćorić, Darga, Briguglio, Schenk, Werner, Andersen and Sames2013). The low relative abundances observed in Morozovelloides here and at Tethys sites (Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; Gebhardt et al. Reference Gebhardt, Ćorić, Darga, Briguglio, Schenk, Werner, Andersen and Sames2013) are therefore likely a result of these subtropical sites being at the ecological limit for the thermophilic Morozovelloides. The biogeographic differences in population dynamics between these two seemingly ecologically similar genera emphasizes the need for spatially replicated ecological sampling.
Stable isotope data, though limited, show that Acarinina and Morozovelloides at Site U1408 had the expected size–δ13C relationship of dinoflagellate symbiont bearers during the MECO (Henehan et al. Reference Henehan, Edgar, Foster, Penman, Hull, Greenop, Anagnostou and Pearson2020). Mixotrophy, or the harboring of photosymbiotic algae, is relatively common in modern planktonic foraminifera (Takagi et al. Reference Takagi, Kimoto, Fujiki, Saito, Schmidt, Kucera and Moriya2019) and has been a key component for shaping spatial and temporal diversity patterns (Ezard et al. Reference Ezard, Aze, Pearson and Purvis2011; Fenton et al. Reference Fenton, Pearson, Jones and Purvis2016b; Hannisdal et al. Reference Hannisdal, Haaga, Reitan, Diego and Liow2017). Despite its continual occurrence throughout geological time in this study, we classify mixotrophy as a specialist, adaptive ecological feeding strategy, as it limits the planktonic foraminifera to a narrow ecological niche (Raia et al. Reference Raia, Carotenuto, Mondanaro, Castiglione, Passaro, Saggese, Melchionna, Serio, Alessio, Silvestro and Fortelius2016; Rolland and Salamin Reference Rolland and Salamin2016). During the middle Eocene, mixotrophic foraminifera likely included Acarinina, Morozovelloides, Globigerinatheka, and Orbulinoides, all of which experienced major global changes in their relative abundance and ecology as a result of transient climate change (Keller Reference Keller1983; Boersma and Silva Reference Boersma and Silva1991; Wade Reference Wade2004; Wade and Pearson Reference Wade and Pearson2008; Wade and Olsson Reference Wade and Olsson2009; Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; Boscolo Galazzo et al. Reference Boscolo Galazzo, Giusberti, Luciani and Thomas2013; Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013), with all (barring several small Acarinina species) becoming extinct before the end of the Eocene (Wade Reference Wade2004; Wade and Pearson Reference Wade and Pearson2008).
Based on the shared ecological strategy of Acarinina and Morozovelloides, one conclusion may be that the reduction of these species at our site was a result of their specialist ecology and changes to their symbiotic relationship (Wade Reference Wade2004; Wade et al. Reference Wade, Al-Sabouni, Hemleben and Kroon2008; Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013), as shown by reduction in test size–δ13C relationship following the MECO at Blake Nose, Site 1051, in the northwest Atlantic (Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013), that occurred as non–symbiont bearing surface layer dwellers continued to thrive (Fig. 4). Yet, despite sharing a similar specialist mixotrophic ecology, Globigerinatheka shows a peak in relative abundance of 33%–34% in the >180 μm size fraction and 11% in the >63 μm size fraction at ~40.40 Ma coincident with peak MECO warming (Fig. 4). In addition, other global records reflect dominance or relative abundance increases of Globigerinatheka through the MECO (Boersma and Premoli Silva Reference Boersma and Silva1986; Boersma et al. Reference Boersma, Silva and Shackleton1987; Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013; Galazzo et al. Reference Galazzo, Thomas, Pagani, Warren, Luciani and Giusberti2014).
Specialist feeding ecologies have been cited as the reason for extinction in deep time of herbivorous sea urchins (Smith and Jeffery Reference Smith and Jeffery1998), herbivorous insects (Labandeira et al. Reference Labandeira, Johnson and Wilf2002), hypercarnivorous canids (Van Valkenburgh Reference Van Valkenburgh2004), and crinoids (Baumiller Reference Baumiller1993). A similar pattern of feeding specialist extinction has also been documented in planktonic foraminifera (Norris Reference Norris1992), but Norris defined a specialist as foraminifera that has limited food sources. The success of Globigerinatheka and persistence of Acarinina and Morozovelloides suggests that specialization is not always entirely detrimental for organisms during transient climatic changes, even with large fluctuations in climate state. The decline in symbiont-bearing planktonic foraminifera across the middle Eocene does suggest the climatic fluctuations pushed these genera closer to their ecological limits (Ezard et al. Reference Ezard, Aze, Pearson and Purvis2011; Edgar et al. Reference Edgar, Bohaty, Gibbs, Sexton, Norris and Wilson2013), which was a process exaggerated at our sites due to the relatively high latitude locality near the species’ biogeographic range limits.
Divergent Response of Size Fraction to the Middle Eocene
Assemblage studies are often conducted at size fractions above >150 μm to avoid juvenile specimens (Al-Sabouni et al. Reference Al-Sabouni, Kucera and Schmidt2007), yet this coarse filter can remove large amounts of diversity and bias studies toward larger individuals, particularly at higher latitudes, where taxa are known to be smaller (Schmidt et al. Reference Schmidt, Thierstein and Bollmann2004). In addition, sampling at a biotically uninformative size fraction can impact inferences on how communities respond to background, transient, and rapid environmental fluctuations. In this study, we found different timings of assemblage responses to middle Eocene climate as a result of size fraction (Figs. 2, 3).
At the relatively high latitude position of our site, water-column heterogeneity was already low due to the general lack of a substantial thermocline at higher latitudes (Rutherford et al. Reference Rutherford, D'Hondt and Prell1999; Al-Sabouni et al. Reference Al-Sabouni, Kucera and Schmidt2007). Background Eocene cooling (Westerhold et al. Reference Westerhold, Marwan, Drury, Liebrand, Agnini, Anagnostou, Barnet, Bohaty, Vleeschouwer, Florindo, Frederichs, Hodell, Holbourn, Kroon, Lauretano, Littler, Lourens, Lyle, Pälike, Röhl, Tian, Wilkens, Wilson and Zachos2020) would have increased water-column stratification, allowing for an increase in relative abundance of genera as a result of widening ecological niches (Whittaker et al. Reference Whittaker, Willis and Field2001; Al-Sabouni et al. Reference Al-Sabouni, Kucera and Schmidt2007). We see the effects of increasing thermal stratification in the larger size fraction (>180 μm), where generic and morphological H S and H GS increase at 42.20 Ma (Fig. 2B,C,E,F). Long-term cooling of the Eocene coupled with thermal stratification and cooling increase before the MECO (Arimoto et al. Reference Arimoto, Nishi, Kuroyanagi, Takashima, Matsui and Ikehara2020; Kearns et al. Reference Kearns, Bohaty, Edgar, Nogué and Ezard2021) results in the removal of larger symbiont-bearing foraminifera (Acarinina and Morozovelloides) and a decline in generic and morphological H S and H GS (Fig. 3). These results imply that amplitude and intensity of environmental change has a major role on how ecosystems respond, possibly larger than the direction of change (Gibbs et al. Reference Gibbs, Bown, Murphy, Sluijs, Edgar, Pälike, Bolton and Zachos2012; Garcia et al. Reference Garcia, Cabeza, Rahbek and Araújo2014; Mayfield et al. Reference Mayfield, Langdon, Doncaster, Dearing, Wang, Velle, Davies and Brooks2021).
In contrast, we do not see any consistent changes in effective diversity at multiple levels of q (Fig. 3) and no substantive differences between pre-MECO and MECO intervals (Fig. 3, Supplementary Table 25). Instead, effective diversity shows significant change in the post-MECO interval (Fig. 3, Supplementary Table 25) with a decrease in morphological, generic, and depth-habitat effective diversity at all levels (q = 0–2) and decreasing assemblage evenness. This trend to less-even communities follows the removal of large symbiont-bearing forms (Fig. 4) and an increase in thermocline dwellers at the expense of mixed-layer species (Fig. 4)
We observe no impact of general Eocene cooling or enhanced pre-MECO cooling on traditional diversity measures in the smaller size fraction compared with the >180 μm size fraction (Fig. 2). In addition, we see no impact of pre-MECO cooling on effective diversity (Fig. 3, Supplementary Table 25). Instead, we observe peaks on morphological and generic richness coinciding with peak MECO warming (Fig. 2A,D). Though the magnitude of warming experienced during the MECO at the sites drilled on IODP Expedition 342 is debated (Arimoto et al. Reference Arimoto, Nishi, Kuroyanagi, Takashima, Matsui and Ikehara2020; Kearns et al. Reference Kearns, Bohaty, Edgar, Nogué and Ezard2021), a global surface ocean temperature increase, alongside the removal of key large symbiont-bearing planktonic foraminifera, may have increased the number of vacant ecological niches, leading to increases in rare, small, microperforate surface-dwelling taxa alongside increases in thermocline dwellers. As a result of the transient nature of the MECO, which lasted ~270–500 kyr (Bohaty and Zachos Reference Bohaty and Zachos2003; Bohaty et al. Reference Bohaty, Zachos, Florindo and Delaney2009; Westerhold and Röhl Reference Westerhold and Röhl2013; Rivero-Cuesta et al. Reference Rivero-Cuesta, Westerhold, Agnini, Dallanave, Wilkens and Alegret2019; Edgar et al. Reference Edgar, Bohaty, Coxall, Bown, Batenburg, Lear and Pearson2020), test size increases in response to increasing warmth, and thus emergence of potentially ecologically optimum conditions, are not observed, unlike at other periods in geological history (Schmidt et al. Reference Schmidt, Renaud and Bollmann2003; Al-Sabouni et al. Reference Al-Sabouni, Kucera and Schmidt2007; Todd et al. Reference Todd, Schmidt, Robinson and Schepper2020). A lack of size response as a result of decreasing thermal stratification across the MECO at this site (Arimoto et al. Reference Arimoto, Nishi, Kuroyanagi, Takashima, Matsui and Ikehara2020; Kearns et al. Reference Kearns, Bohaty, Edgar, Nogué and Ezard2021) was a driver of globally small test sizes in the middle Eocene (Schmidt et al. Reference Schmidt, Thierstein and Bollmann2004) and was responsible for the emergence of a latitudinal size gradient at ~42 Ma that persists today (Schmidt et al. Reference Schmidt, Thierstein and Bollmann2004). The brevity of the MECO also meant increasing diversity in the >63 μm size fraction was short-lived and was followed by a dramatic decline in generic and morphological richness (Fig. 3A,D) and effective diversity (q <1; Fig. 3). Post-MECO cooling had little other effect on effective diversity, except a reduction in morphological diversity (q = 1–2) as a result of the decline in morphologically distinct Chiloguembelina and Jenkinsina.
Insights into Paleoceanographic Changes across the MECO from “Rare” Taxa
In our record, we have numerous genera with low occurrence (~8 genera) that could be considered rare and cause large fluctuations in generic and morphogroup effective diversity (q = 0) at both size fractions (Figs. 2, 3). Although rarity in itself is potentially an important measure (e.g., acting as canaries for early warning signals; Doncaster et al. Reference Doncaster, Chávez, Viguier, Wang, Zhang, Dong, Dearing, Langdon and Dyke2016), being rare is common, with the majority of taxa represented by only a few individuals (Gaston Reference Gaston2008). Microperforate biserial and triserial taxa such as Chiloguembelina and Jenkinsina are rare in many records, as they occur in highest abundance in the infrequently studied >63 μm size fraction despite being omnipresent throughout the Cenozoic (Li and Radford Reference Li and Radford1991), with approximately 20 species occurring in the Eocene alone (Huber et al. Reference Huber, Olsson, Pearson, Pearson, Olsson, Huber, Hemleben and Berggren2006). In addition, these taxa have sporadic geographic and biostratigraphic records (Kroon and Nederbragt Reference Kroon and Nederbragt1990; Darling et al. Reference Darling, Thomas, Kasemann, Seears, Smart and Wade2009), often increasing to noticeable abundances during periods of environmental stress such as ocean acidification events (Nederbragt et al. Reference Nederbragt, Erlich, Fouke and Ganssen1998; Coccioni et al. Reference Coccioni, Luciani and Marsili2006), periods of background climatic instability (Kroon and Nederbragt Reference Kroon and Nederbragt1990; Li and Radford Reference Li and Radford1991; Luciani et al. Reference Luciani, Giusberti, Agnini, Backman, Fornaciari and Rio2007, Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010; D'Haenens et al. Reference D'Haenens, Bornemann, Roose, Claeys and Speijer2012), and the aftermath of mass extinction events (Keller Reference Keller1993; Luciani Reference Luciani1997, Reference Luciani2002; Keller et al. Reference Keller, Adatte, Stinnesbeck, Luciani, Karoui-Yaakoub and Zaghbib-Turki2002). The lack of changes at the >63 μm size fraction for most diversity measures compared with the >180 μm size fraction supports these results that smaller taxa are more resilient to both background (i.e., Eocene cooling) and transient perturbations (MECO) compared with large taxa and thus deserving of further study.
In our record, both Chiloguembelina and Jenkinsina are only substantive components of assemblages in the >63 μm size fraction (Fig. 4), not at >180 μm. Two noticeable peaks occur in Chiloguembelina (43.24% at 39.85 Ma and 52.27% at 41.45 Ma) and Jenkinsina (30.42% at 40.41 Ma and 20.05% at 41.45 Ma), coinciding with paleoceanographic instability across the MECO, interrupting otherwise relatively low relative abundances (Chiloguembelina: <~20%, Jenkinsina: <~1%; Fig. 4). Similar peaks in abundance of Chiloguembelina and Jenkinsina have been observed in the Tethys Ocean (Alano section; Luciani et al. Reference Luciani, Giusberti, Agnini, Fornaciari, Rio, Spofforth and Pälike2010) and at other high-latitude sites (Li and Radford Reference Li and Radford1991) and associated with upwelling or low-oxygen conditions. There is no evidence for either at our study site (Arimoto et al. Reference Arimoto, Nishi, Kuroyanagi, Takashima, Matsui and Ikehara2020; Kearns et al. Reference Kearns, Bohaty, Edgar, Nogué and Ezard2021), but these peaks in abundance support arguments that these taxa thrive during transient climatic events.
One hypothesis as to why rare taxa flourish during environmental perturbations is that they replace superior predators or dominant taxa that are lost in order to maintain ecological function (Walker et al. Reference Walker, Kinzig and Langridge1999). Both Chiloguembelina relative abundance peaks occur synchronously with troughs in Acarinina relative abundance (Fig. 4). As all three genera (Chiloguembelina, Morozovelloides, and Acarinina) are surface dwelling, this rise in abundance may be a response to the sparsely occupied ecological niche(s) left by the removal of a large proportion of Acarinina and Morozovelloides. This synchronicity may reflect changes in surface-water productivity, as Chiloguembelina is an opportunistic eutrophic genus (Luciani et al. Reference Luciani, D'Onofrio, Filippi and Moretti2020). Coincidental changes in relative abundance also occur in the thermocline and subthermocline, where Jenkinsina peaks at the same time as Subbotina (Fig. 4). While no taxa are removed from the thermocline during the MECO, Subbotina has a wide and plastic ecological niche (Kearns et al. Reference Kearns, Bohaty, Edgar, Nogué and Ezard2021), and it may be that Jenkinsina prospered for a short interval due to temporary availability of the thermocline ecological niche. Stable isotope studies of Chiloguembelina and Jenkinsina across the middle Eocene would be one route to rigorously testing this hypothesis. Our observations do suggest, however, that these rara taxa are useful when looking at paleoceanographic changes to identify periods of environmental instability and therefore should be measured instead of being dismissed for their small size and sporadic geographic and biostratigraphic records (Kroon and Nederbragt Reference Kroon and Nederbragt1990; Darling et al. Reference Darling, Thomas, Kasemann, Seears, Smart and Wade2009).
Conclusion
Our analysis of planktonic foraminiferal assemblages within the middle Eocene at sites on the Newfoundland margin demonstrate that complex diversity dynamics follow transient environmental changes. We show that transient events are not necessarily terminal for specialist taxa but can push these taxa to their ecological limits, which potentially influences their abundance and community composition for millions of years after a transient perturbation. We argue that rather than complicating our understanding of planktonic foraminifera responses in the middle Eocene, measuring at two size fractions illuminates size dynamics more fully to enhance our understanding of paleooceanographic drivers of biotic turnover. The impact of rare taxa is not well described by its standing percentage in a community. However, by documenting the smaller size fraction, we were able to record smaller and more microperforate taxa, not normally measured, to show how rare taxa can provide equivalent ecological function to those taxa that are lost and inform our understanding of environmental perturbations and how communities to persist through climate fluctuations.
Acknowledgments
This research used samples provided by the IODP. We would like to thank A. Brombacher and M. Holmström for taxonomic guidance. We would also like to thank M. Yasuhara and one anonymous reviewer for their helpful comments when reviewing this article. This work was supported by the Natural Environment Research Council awards NE/L002531/1 and NE/P019269/1, as well as the Institute of Life Sciences (IfLS) at the University of Southampton. The authors declare no competing interests.
Data Availability Statement
Data available from the Dryad digital repository: https://doi.org/10.5061/dryad.cfxpnvx81.
Open access
