Small but mighty: how overlooked small species maintain community structure through middle Eocene climate change

Abstract Understanding current and future biodiversity responses to changing climate is pivotal as anthropogenic climate change continues. This understanding is complicated by the multitude of available metrics to quantify dynamics and by biased sampling protocols. Here, we investigate the impact of sampling protocol strategies using a data-rich fossil record to calculate effective diversity using Hill numbers for the first time on Paleogene planktonic foraminifera. We sample 22,830 individual tests, in two different size classes, across a 7 Myr time slice of the middle Eocene featuring a major transient warming event, the middle Eocene climatic optimum (MECO; ~40 Ma), at study sites in the midlatitude North Atlantic. Using generalized additive models, we investigate community responses to climatic fluctuations. After correcting for any effects of fossil fragmentation, we show a peak in generic diversity in the early and middle stages of the MECO as well as divergent trajectories between the typical size-selected community (>180 μm) and a broader assemblage, including smaller genera (>63 μm). Assemblages featuring smaller genera are more resilient to the climatic fluctuations of the MECO than those assemblages that feature only larger genera, maintaining their community structure at the reference Hill numbers for Shannon's and Simpson's indices. These results raise fundamental questions about how communities respond to climate excursions. In addition, our results emphasize the need to design studies with the aim of collecting the most inclusive data possible to allow detection of community changes and determine which species are likely to dominate future environments.


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. 1994;Purvis and Hector 2000;Mace et al. 2012;Hohenegger 2014), 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 2000) understandable to both researchers and the general public (Purvis and Hector 2000;Baum 2009;Chiarucci et al. 2011;Reydon 2019), 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 1942) and share phenotypic and ecological traits . 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. 2014). 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. 2014).
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. 2017), and are also likely to be influential components of ecosystems (Lennon et al. 2004;Gaston 2008;Hannisdal et al. 2017). Thus, it is hypothesized that common taxa contribute more to assemblage diversity than rare taxa (Lennon et al. 2004;Gaston 2008). 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 2008), 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 2007;Gaston 2008;Hannisdal et al. 2017), 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. 2016) and insights into paleoceanographic change (Ishino and Suto 2020) and acting as the focus of conservation efforts (Gaston 2008). However, taxa vary spatially, directly influencing ecological dynamics and diversity patterns within each generation at any single location only (Patzkowsky and Holland 2007). Consequently, what is rare in one sample or area may be common in another (Colwell 2009). 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. 2011). On a theoretical level, far less is known about the functional role of rare taxa in their ecosystems (Lyons et al. 2005), meaning they are easier to dismiss as unimportant and therefore to ignore (Chao et al. 2014b).
To be common, rare, or absent, is a relative measure (Preston 1948), 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 2000;Colwell 2009;Reich et al. 2012) 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 2010). Such methods have been applied: effective numbers, or Hill numbers (Hill 1973), 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. 2014b), spiders (Chao et al. 2014b(Chao et al. , 2020, corals (Chao et al. 2020), and bacteria (Kang et al. 2016). In the fossil record, assemblages are at the mercy of time and preservation (Jackson and Blois 2015). 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. 2021), deep-sea ostracods (Yasuhara et al. 2016), and pelagic planktonic foraminifera (Yasuhara et al. 2020), as well as Paleozoic marine radiations (Rasmussen et al. 2019).
Conceptually, Hill numbers are the effective number of equally abundant taxa required to give the same diversity presented in the sample (Hill 1973;Jost 2010a;Chao et al. 2014bChao et al. , 2020. 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 1973). 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 1973;Chiu and Chao 2014). 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. 2019;Miller et al. 2020;Heritier-Robbins et al. 2021) or, as in the present study, geologically transient 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. climatic events. The commonality of units at all levels of q means that inferences can be made regarding magnitudes of change (Jost 2007(Jost , 2010aChao et al. 2014b) and sample and locality differences (Hill 1973;Chao et al. 2014a) 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. 2010), functional (Chiu and Chao 2014), and taxonomic (Chao et al. 2014a) diversity, with straightforward bootstrapping techniques to quantify how high proportions of singletons increase uncertainty (Chao et al. 2014a).
To be meaningful, Hill number calculations require sufficient and careful sampling protocols. Here we calculate Hill numbers for a deeptime 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 2003;Bohaty et al. 2009;Rivero-Cuesta et al. 2019;Edgar et al. 2020). The requirement for large samples of individuals means fossilized planktonic foraminifera are an ideal candidate for Hill numbers (Yasuhara et al. 2020) 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 2001;Kucera 2007;Brummer and Kučera 2022), which upon death are deposited on the seafloor in vast quantities, producing ∼2 Gt of calcite per year (Schiebel and Hemleben 2008). Deposition of planktonic foraminifera has occurred nearly continuously since their evolution ∼200 Ma during the Jurassic period (Fraass et al. 2015), 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 Fraass et al. 2015;Fenton et al. 2016a;Yasuhara et al. 2017) with a highly temporally and spatially resolved record (Fenton et al. 2021), 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. 2018Westerhold et al. , 2020 through to the "icehouse" of the Oligocene that started at the Eocene/Oligocene transition (EOT, 34 Ma;Westerhold et al. 2020;Hutchinson et al. 2021) with the establishment of continent-wide Antarctic glaciation (Zachos et al. 1996;Coxall et al. 2005). 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 2003;Bohaty et al. 2009;Rivero-Cuesta et al. 2019). During the MECO, there was a transient ∼3°C-6°C rise in surface and deep-water temperatures (Bohaty and Zachos 2003;Bohaty et al. 2009;Bijl et al. 2010;Galazzo et al. 2014;Cramwinckel et al. 2019;Henehan et al. 2020), reduced surface ocean pH (Henehan et al. 2020), and a shoaling of the calcium carbonate compensation depth (CCD; Bohaty and Zachos 2003;Bohaty et al. 2009). The MECO is terminated by a rapid return to pre-MECO conditions (Bohaty et al. 2009) 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 1971;Boersma and Premoli Silva 1986;Boersma and Silva 1991;Keller et al. 1992;Wade 2004;Sexton et al. 2006;Wade and Pearson 2008;Luciani et al. 2010;Ezard et al. 2011;Galazzo et al. 2014;Fenton et al. 2016a). Middle Eocene biotic changes in planktonic foraminifera include: (1) the progressive extinction of surface-dwelling, symbiont-bearing taxa (Boersma and Premoli Silva 1986;Boersma and Silva 1991;Keller et al. 1992;Wade 2004;Wade and Pearson 2008); (2) a reduction in test size (Schmidt et al. 2004;Wade and Pearson 2008;Wade and Olsson 2009); (3) development of latitudinal size (Schmidt et al. 2004) and diversity (Fenton et al. 2016a) gradients alongside major assemblage fluctuations (Steineck 1971;Keller 1983;Boersma and Premoli Silva 1986;Boersma et al. 1987;Hallock et al. 1991;Keller et al. 1992;Sexton et al. 2006;Luciani et al. 2010;Galazzo et al. 2014); and (4) changes in ecology, for example, loss or inhibition of algal photosymbionts from hosting taxa Edgar et al. 2013) and shallowing depth habitat of Hantkenina (Coxall et al. 2000). 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. 2008;Wade and Pearson 2008;Pearson and Wade 2015). The MECO resulted in a global crisis for muricate taxa (Luciani et al. 2010;Edgar et al. 2013), varying symbiotic taxon responses (Luciani et al. 2010;Edgar et al. 2013;Gebhardt et al. 2013;Arimoto et al. 2020;Kearns et al. 2021), and increased abundance of ecologically flexible (Galazzo et al. 2015;Kearns et al. 2021) and small opportunistic taxa (Luciani et al. 2010). 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.
Sample ages from Sites U1408 and U1410 were calculated based on age-depth models constructed using available biostratigraphy and magnetostratigraphy (Norris et al. 2014). The 2012 geological timescale was then used for age calibrations for the middle Eocene geomagnetic reversals (GTS2012; Gradstein et al. 2012). Sample ages for Site U1406 are based upon shipboard biostratigraphic and magnetostratigraphic data (Norris et al. 2014;Van Peer 2017). 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. 2021). 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. 2005) to avoid sampling juvenile specimens and to enable species-level identification (Al-Sabouni et al. 2007. 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. 2007). 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. 1986).

Diversity Analysis
A sample of 300 individuals is considered sufficient to estimate diversity in foraminiferal assemblages (Al-Sabouni et al. 2007), despite the potential for missing rare specimens due to low abundances (Jost 2010b). 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 Wade et al. 2018).
To understand diversity changes further, we then classified all genera into morphogroups (Supplementary Tables 5, 6) adapted from previous classifications ) 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. 2006;Wade et al. 2018). 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 1973): 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 2012; 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 2012; 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 (2015), in the Supplementary Material. Effective diversity was calculated at the default 0.1 intervals for q between 0 and 2 (Chao and Jost 2015; 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 1971;Malmgren 1987;Nguyen et al. 2009). 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. 2009Nguyen et al. , 2011Nguyen and Speijer 2014) and is influenced by the microenvironment of the individual, which influences test chemistry and causes interspecific differences in dissolution susceptibility (Berger 1970;Nguyen et al. 2011;Petro et al. 2018). To account for this variability, we use an accepted fragmentation proxy to estimate the dissolution levels (Le and Shackleton 1992), using the proportion of planktonic foraminiferal test fragments (Frag) and whole specimens: (2) We classify a fragment as anything <75% of a whole specimen (which is more conservative than the <50% previously used; Malmgren 1987). Foraminifera have a tendency to break into multiple pieces; therefore, the percentage of fragments in a sample varies nonlinearly with dissolution (Le and Shackleton 1992). 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 1992; Leon-Rodriguez and Dickens 2010). We use a baseline of 20% fragmentation to indicate normal levels of fragmentation and dissolution (Pfuhl and Shackleton 2004). 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 2017) 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 crossvalidation 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 2011). 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).
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. 2022) was applied, due to unequal observations in each interval (Zar 2010), to identify intervals which were significantly different from each other.

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 2004), 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 backtransformed 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 , implying that the second-ranked models had "essentially no" support (Burnham and Anderson 2002).
Following model selection, we used the bestfitting model (Tables 1, 2) to predict diversity values across our study interval at a mean fragmentation of 10% based on the mean of our   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.
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).

Hill Numbers and Relative Abundance Fluctuations
Genera.-When relative abundance is not considered (q < 1), pre-MECO (>41.09 Ma),    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).
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 depthhabitat 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).
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  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. 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 1971;Thunell 1976), biasing climatic and biotic interpretations of these assemblages (Berger 1973). Dissolution can be morphologically selective (Berger 1970;Boltovskoy and Totah 1992;Petrizzo et al. 2008;Nguyen et al. 2009, 2011, but see Petro et al. 2018) with speciesspecific tendencies (Nguyen et al. 2011;but see Berger 1970;Malmgren 1987). Across the MECO, extensive dissolution as a result of shoaling CCD has been recorded in the Pacific, Indian, and Atlantic Oceans (Lyle et al. 2005;Bohaty et al. 2009;Pälike et al. 2012). Despite our study sites sitting well above the known late Paleogene CCD (Norris et al. 2014), 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 2004). 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. 2013;Savian et al. 2014), there was no observable drop in overall planktonic foraminifera abundances, which would indicate dissolution impacted assemblages (Malmgren 1987). We do detect a statistically significant effect of dissolution in 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). 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 openocean 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. 2013); 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. 2010). 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. 2010;Fig. 4), unlike the lowerlatitude Blake Nose (∼25°N), where Acarinina recovers in both abundance and test size (Edgar et al. 2013). 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 Luciani et al. 2010;Edgar et al. 2013), although in the Tethys, Morozovelloides is scarcely abundant throughout the middle Eocene (Luciani et al. 2010;Gebhardt et al. 2013). The low relative abundances observed in Morozovelloides here and at Tethys sites (Luciani et al. 2010;Gebhardt et al. 2013) 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-δ 13 C relationship of dinoflagellate symbiont bearers during the MECO (Henehan et al. 2020). Mixotrophy, or the harboring of photosymbiotic algae, is relatively common in modern planktonic foraminifera (Takagi et al. 2019) and has been a key component for shaping spatial and temporal diversity patterns Fenton et al. 2016b;Hannisdal et al. 2017). 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. 2016;Rolland and Salamin 2016). 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 1983;Boersma and Silva 1991;Wade 2004;Wade and Pearson 2008;Wade and Olsson 2009;Luciani et al. 2010;Boscolo Galazzo et al. 2013;Edgar et al. 2013), with all (barring several small Acarinina species) becoming extinct before the end of the Eocene (Wade 2004;Wade and Pearson 2008).
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 2004;Edgar et al. 2013), as shown by reduction in test size-δ 13 C relationship following the MECO at Blake Nose, Site 1051, in the northwest Atlantic (Edgar et al. 2013), 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 1986;Boersma et al. 1987;Edgar et al. 2013;Galazzo et al. 2014).
Specialist feeding ecologies have been cited as the reason for extinction in deep time of herbivorous sea urchins (Smith and Jeffery 1998), herbivorous insects (Labandeira et al. 2002), hypercarnivorous canids (Van Valkenburgh 2004), and crinoids (Baumiller 1993). A similar pattern of feeding specialist extinction has also been documented in planktonic foraminifera (Norris 1992), 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 Edgar et al. 2013), 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. 2007), 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. 2004). 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. 1999;Al-Sabouni et al. 2007). Background Eocene cooling (Westerhold et al. 2020) 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. 2001;Al-Sabouni et al. 2007). 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. 2020;Kearns et al. 2021) 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. 2012;Garcia et al. 2014;Mayfield et al. 2021).
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. 2020;Kearns et al. 2021), 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 2003;Bohaty et al. 2009;Westerhold and Röhl 2013;Rivero-Cuesta et al. 2019;Edgar et al. 2020), 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. 2003;Al-Sabouni et al. 2007;Todd et al. 2020). A lack of size response as a result of decreasing thermal stratification across the MECO at this site (Arimoto et al. 2020;Kearns et al. 2021) was a driver of globally small test sizes in the middle Eocene (Schmidt et al. 2004) and was responsible for the emergence of a latitudinal size gradient at ∼42 Ma that persists today (Schmidt et al. 2004). 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. 2016), being rare is common, with the majority of taxa represented by only a few individuals (Gaston 2008). 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 1991), with approximately 20 species occurring in the Eocene alone . In addition, these taxa have sporadic geographic and biostratigraphic records (Kroon and Nederbragt 1990;Darling et al. 2009), often increasing to noticeable abundances during periods of environmental stress such as ocean acidification events (Nederbragt et al. 1998;Coccioni et al. 2006), periods of background climatic instability (Kroon and Nederbragt 1990;Li and Radford 1991;Luciani et al. 2007Luciani et al. , 2010D'Haenens et al. 2012), and the aftermath of mass extinction events (Keller 1993;Luciani 1997Luciani , 2002Keller et al. 2002). 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. 2010) and at other high-latitude sites (Li and Radford 1991) and associated with upwelling or low-oxygen conditions. There is no evidence for either at our study site (Arimoto et al. 2020;Kearns et al. 2021), 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. 1999). 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. 2020). 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. 2021), 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 1990;Darling et al. 2009).

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.