Searching for Function: Reconstructing Adaptive Niche Changes Using Geochemical and Morphological Data in Planktonic Foraminifera

Dead species remain dead. The diversity record of life is littered with examples of declines and radiations, yet no species has ever re-evolved following its true extinction. In contrast, functional traits can transcend diversity declines, often develop iteratively and are taxon-free allowing application across taxa, environments and time. Planktonic foraminifera have an unrivaled, near continuous fossil record for the past 200 million years making them a perfect test organism to understand trait changes through time, but the functional role of morphology in determining habitat occupation has been questioned. Here, we use single specimen stable isotopes to reconstruct the water depth habitat of individual planktonic foraminifera in the genus Subbotina alongside morphological measurements of the tests to understand trait changes through the Middle Eocene Climatic Optimum [MECO: ∼40 Myr ago (mega annum, Ma)]. The MECO is a geologically transient global warming interval that marks the beginning of widespread biotic reorganizations in marine organisms spanning a size spectrum from diatoms to whales. In contrast to other planktonic foraminiferal genera, the subbotinids flourished through this interval despite multiple climatic perturbations superimposed on a changing background climate. Through coupled trait and geochemical analysis, we show that Subbotina survival through this climatically dynamic interval was aided by trait plasticity and a wider ecological niche than previously thought for a subthermocline dwelling genus supporting a generalist life strategy. We also show how individually resolved oxygen isotopes can track shifts in depth occupancy through climatic upheaval. During and following the MECO, temperature changes were substantial in the thermocline and subthermocline in comparison to the muted responses of the surface ocean. In our post-MECO samples, we observe restoration of planktonic foraminifera depth stratification. Despite these changing temperatures and occupied depths, we do not detect a contemporaneous morphological response implying that readily available traits such as test size and shape do not have a clear functional role in this generalist genus. Modern imaging measurement technologies offer a promising route to gather more informative morphological traits for functional analysis, rather than the traditional candidates that are most easily measured.

Dead species remain dead. The diversity record of life is littered with examples of declines and radiations, yet no species has ever re-evolved following its true extinction. In contrast, functional traits can transcend diversity declines, often develop iteratively and are taxon-free allowing application across taxa, environments and time. Planktonic foraminifera have an unrivaled, near continuous fossil record for the past 200 million years making them a perfect test organism to understand trait changes through time, but the functional role of morphology in determining habitat occupation has been questioned. Here, we use single specimen stable isotopes to reconstruct the water depth habitat of individual planktonic foraminifera in the genus Subbotina alongside morphological measurements of the tests to understand trait changes through the Middle Eocene Climatic Optimum [MECO: ∼40 Myr ago (mega annum, Ma)]. The MECO is a geologically transient global warming interval that marks the beginning of widespread biotic reorganizations in marine organisms spanning a size spectrum from diatoms to whales. In contrast to other planktonic foraminiferal genera, the subbotinids flourished through this interval despite multiple climatic perturbations superimposed on a changing background climate. Through coupled trait and geochemical analysis, we show that Subbotina survival through this climatically dynamic interval was aided by trait plasticity and a wider ecological niche than previously thought for a subthermocline dwelling genus supporting a generalist life strategy. We also show how individually resolved oxygen isotopes can track shifts in depth occupancy through climatic upheaval. During and following the MECO, temperature changes were substantial in the thermocline and subthermocline in comparison to the muted responses of the surface ocean. In our post-MECO samples, we observe restoration of planktonic foraminifera depth stratification. Despite these changing temperatures and occupied depths, we do not detect a contemporaneous morphological response implying that readily available

INTRODUCTION
Existence in an ecosystem implies importance. Abundance estimates and extinction rates are key measures used to monitor efforts to conserve species but counting alone cannot accurately demonstrate the health of an ecosystem and the species within it (Akçakaya et al., 2020). To be present in an ecosystem an organism must have a role, and thus a function (Jax, 2005), which is fluid through time, space and biota (Akçakaya et al., 2020). The key is therefore not to merely note an organism's presence but to understand the traits that allows it to function within an ecosystem such as pollinator body and face hairiness which aids pollination (Stavert et al., 2016) or wing length which influences foraging distance (Brousseau et al., 2018). Functional loss through species extinction can be undetectable, if another species has the same or a similar function, or non-linear and saturating (Cardinale et al., 2012), with the smallest of functions having large impacts on ecosystem health (Akçakaya et al., 2020). But function is not yet used effectively in conservation. For example, functional rarity is highest for species not identified at risk on the IUCN red list in global coral reefs (Grenié et al., 2018) and simulations based on Californian bee populations show the most functionally efficient and important bees (from an analysis of 12 species/genera) are most at risk of extinction (Larsen et al., 2005). Therefore, identifying and understanding function is of paramount importance.
Simply described, a function typically describes some form of ecological process (Farnsworth et al., 2017) yet applying this definition in ecology has led to much ambiguity (Jax, 2005;Farnsworth et al., 2017). In modern ecosystems, assigning function to a trait can be done through observation and experiments with direct reference to human requirements and usefulness (Mace et al., 2014). In situ observations are one dimensional, however, and if the past is the key to predicting the future (Tierney et al., 2020), then we need to test modes of inference of how organisms lived thousands to millions of years ago and we need to be able to integrate modern and fossil functional indices. In the fossil record, direct observations of physiology and ecological function are typically not possible, so we are often left instead to infer an organism's functional role in its community from the preserved morphological traits. Such traits may therefore be a viable currency to "bind the past and present together" (Eronen et al., 2010). Here, we relate morphological traits in fossilized planktonic foraminifera to inferred depth habitats, derived from geochemical measurements, through climatic change during the middle Eocene (∼48 to 38 Ma).
We define a trait as any morphological, physiological, phenological or behavioral feature measurable at the individual level (Violle et al., 2007). Traits are the avenue through which an organism interacts with its environment (biotic and abiotic) (Lacourse, 2009;Oksanen et al., 2019) and determine whether an organism survives or meets its demise in a changing environment (McGuire and Lauer, 2020). Thus, traits capture the most important aspects of the environment, and consequently it is the trait and not taxonomic identity that is crucial (Lacourse, 2009). For a trait to be classed as functional, it must, directly or indirectly, impact individual performance and fitness of species (McGill et al., 2006).
Taxon-free traits are a specific class of traits that transcend taxonomic classification and provide a commonality to allow comparisons across communities in different climatic and geographical settings (McGill et al., 2006). Often these traits are phenotypic such as dental morphology (e.g., Renaud et al., 2005;Žliobaitë et al., 2016;Oksanen et al., 2019;McGuire and Lauer, 2020), overall body geometry (Bregman et al., 2014;Pimiento et al., 2017;Macumber et al., 2020;Antczak-Orlewska et al., 2021;Di Martino and Liow, 2021) and, in our study system planktonic foraminifera, test morphology (Baumfalk et al., 1987;Huber et al., 1997;Renaud and Schmidt, 2003;Schmidt et al., 2004Schmidt et al., , 2006Payne et al., 2012;Rego et al., 2012;Weinkauf et al., 2014Weinkauf et al., , 2019Weiner et al., 2015;Brombacher et al., 2017;Kucera et al., 2017). Traits can also include ecological traits like habitat (Bregman et al., 2014;Pimiento et al., 2017) and feeding behavior (Bregman et al., 2014). Despite morphological traits being collected in abundance (Parr et al., 2012), trait based research on fauna is reduced compared to flora (Lavorel and Garnier, 2002;Lacourse, 2009;Steinthorsdottir et al., 2016;Fried et al., 2019;Birks, 2020), which has resulted in a plethora of plant-based data bases [e.g., TRY (Kattge et al., 2020); LEDA (Kleyer et al., 2008) and BIEN (Enquist et al., 2016)]. One reason for this discrepancy, is that faunal morphological traits were collected before the birth of trait-based ecology and are therefore not located in easy-toaccess databases and require advancements and applications of text mining tools (Parr et al., 2012). Another fundamental issue is that plants make up 81% of global biomass compared to the 0.73% made up by protists (including planktonic foraminifera) and 0.37% made up by animals ( Bar-On et al., 2018), meaning flora is more accessible for research. One promising approach then is to apply trait-based methods in the geological record to understand faunal responses to environmental change.
We have the clearest indication of the impact of extinction on species richness/taxon counts in the deep time fossil record, but assignment of a functional trait in extinct species implies some knowledge of the environment (Violle et al., 2007), as well as an observation of how the trait reacts to climatic fluctuations (Eronen et al., 2010). This correlative relationship allows us to infer the (biotic and abiotic) environment from the traits observed, or vice versa if we wish (Eronen et al., 2010). If function only exists within the context of the broader community, then this inference is particularly challenging and prone to inconsistent extrapolation. Ideally, we should measure traits with an inferred ecological function alongside an environmental indicator that indicates habitat, and then seek to match the two signals.
Our goal in this study is to understand trait changes within the context of a broader community undergoing large-scale abiotically induced change. We present new morphological and geochemical data on planktonic foraminifera. Planktonic foraminifera are extant, holoplanktonic, single celled organisms that build calcite shells (tests) that, upon death, rain down in abundance to the seafloor contributing significantly to deepsea biogenic carbonate (Vincent and Berger, 1981). Foraminifera tests provide a near continuous, spatially and temporally highresolution fossil record dating back ∼200 Myrs to the Jurassic (Fraass et al., 2015). Additionally, planktonic foraminifera are currently represented by ∼50 extant species (Kucera, 2007) meaning specimens can be cultured in the lab (e.g., Bé et al., 1981Bé et al., , 1982Bijma et al., 1990;Spero and Lea, 1993;Henehan et al., 2017;Burke et al., 2018;Fehrenbacher et al., 2018;Holland et al., 2020), yielding an ever-growing knowledge of micro and macroscale influences on foraminifera life-history. In addition, our understanding of what are functional traits in foraminifera is increasing. Some traits have obvious function such as spines for feeding (Hemleben et al., 1991) or the presence of symbionts (Bé et al., 1977(Bé et al., , 1981Takagi et al., 2019). Recent studies, building on previous hypotheses, have shown that pores on foraminifera tests may be functionally linked to gas exchange (Bé, 1968;Baumfalk et al., 1987;Burke et al., 2018) but are only measured reliably from the inside of the test requiring scanning electron microscope (SEM) images or micro-CT scanning (Constandache et al., 2013;Burke et al., 2018). In contrast, traits such as test shape and size must have a function, but what those functions are remain highly debated (Renaud and Schmidt, 2003;Caromel et al., 2014;Burke and Hull, 2017).
A high-resolution fossil record, advances in morphological measurements and increasing knowledge of function mean planktonic foraminifera are an ideal candidate to investigate trait-based responses to environmental change in the geological record. Furthermore, geochemical methods are advancing at an astounding rate making ecological inferences such as life history, metabolic rate, gene flow and geolocations accessible in deep time (Trueman et al., 2016). In planktonic foraminifera these advances have allowed geochemical measurements to be taken at the level of the individual, enhancing our knowledge of paleoclimate (Thirumalai et al., 2013;Schmitt et al., 2019;Glaubke et al., 2021) and the impact of individual planktonic foraminifera ecology on geochemical signatures (Eggins et al., 2003(Eggins et al., , 2004Friedrich et al., 2012;Fehrenbacher et al., 2018;Groeneveld et al., 2019;Weinkauf et al., 2020). We leverage these advances in analytical techniques alongside the exemplary evolutionary record of planktonic foraminifera to investigate trait, organismal and community responses to climatic change on geological timescales. Here, we focus on planktonic foraminifera trait changes across a transient warming event known as the Middle Eocene Climatic Optimum (MECO).

Material and Sample Preparation
Each sample was dried at 40 • C for 4 days and then soaked in sodium hexametaphosphate for a minimum of 4 days on a shaker table to disaggregate the sediment. Samples were then washed over a 38 µm sieve using deionized water until the water ran clear and then dried overnight in a 40 • C oven overnight before being transferred to vials. Each sample was subsequently dry sieved to allow picking of individual specimens under the light microscope. For multi-specimen analysis (Section "Stable isotope analysis") foraminifera were picked from a narrow sieve size fraction of 250-315 µm to avoid the effects of foraminifera size on geochemical analysis. For the individual foraminifera analysis (Section "Individual foraminifera analysis") we picked from the > 180 µm size fraction to capture the widest range of morphological variation from the adult population.

Stable Isotope Analysis
Planktonic foraminifera occupy different depth niches within the water column creating an ecologically stratified community. As a result of isotopic fractionation in the water column, the stable isotopic signature of planktonic foraminifera tests reflects the depth habitat in which they live. The fractionation of oxygen into foraminifera calcite is temperature dependent, resulting in increasing foraminifera test δ 18 O values with increasing depth paralleling the trend of decreasing temperature with depth in the oceans (Fairbanks et al., 1980(Fairbanks et al., , 1982. The opposite depth dependent trend is seen in carbon isotopes with foraminifera calcite δ 13 C decreasing with depth. This relationship is due to the preferential uptake of 12 C during photosynthesis at shallow water depths and the export of particulate organic carbon from the upper water column. Remineralization of particulate organic carbon at depth then releases isotopically light carbon back into the dissolved organic carbon (DIC) pool. Thus, surface ocean symbiont-hosting foraminifera have relatively low δ 18 O and high δ 13 C values relative to nonsymbiont subthermocline dwellers with thermocline dwelling sitting somewhere in-between (Figure 1). These relationships do assume isotopic equilibrium between the foraminifera test and seawater. However, this equilibrium is offset by a number of physical parameters such as salinity, carbonate ion concentration (Urey, 1947;Urey et al., 1951;Epstein et al., 1953;Epstein and Lowenstam, 1954;Spero et al., 1997;Pearson, 2012) as well as foraminifera biology and ecology often referred to as "vital effects" (Erez, 1978;Spero and Williams, 1989;Spero et al., 1991;Bemis et al., 1998;Friedrich et al., 2012;Birch et al., 2013). These factors, especially foraminifera biology (Edgar et al., 2017), tend to impact test δ 13 C more than δ 18 O and can be minimized, for example, by picking foraminifera within a narrow size fraction. Nevertheless, the broad patterns of depth ranking using δ 18 O and δ 13 C (Figure 1) remains true and can be used to reconstruct foraminifera depth habitats through geological time  Pearson (1998). (Pearson, 1998;Spero, 1998;Pearson et al., 2001;Sexton et al., 2006a,b;Coxall et al., 2007;Birch et al., 2013). To investigate water column structure across our study interval, planktonic foraminifera from genera representing three depth ecologies were picked: Globigerinatheka (mixed layer), Subbotina (thermocline) and Catapsydrax (sub-thermocline) (Figure 1).
We use the symbiont bearing, deep mixed layer dwelling Globigerinatheka (Sexton et al., 2006b;Edgar et al., 2013) to represent the mixed layer (Figure 1) rather than other commonly used shallower mixed layer inhabitants Acarinina or Morozovelloides. Morozovelloides which were not abundant enough in our samples for statistically robust isotope or morphological analyses, whilst Acarinina were abundant but showed signs of reworking (stained and heavily fragmented) and recrystallization. We made every effort to pick individuals with no visible signs of either gametogenic calcite or recrystallized wall textures and were mindful of the caveats during interpretation. To represent the thermocline and subthermocline we used asymbiotic Subbotina and Catapsydrax, respectively (Figure 1).
Approximately 20 individuals from each of Globigerinatheka spp., Subbotina spp. (S. utilisindex, S. eoceana, S. projecta, and S. linaperta) and Catapsydrax unicavus were picked from each of the 6 samples following the taxonomy of Pearson et al. (2006). The 18 subsamples, each of 20 individuals, were then crushed, homogenized and weighed into vials. Each subsample weighed between 50 and 60 µg (Supplementary Table 2) and was then cleaned by ultrasonification in ethanol for 3-5 s, rinsed in deionized water and then placed in a 40 • C oven for 1-2 h to dry. Coupled δ 18 O and δ 13 C measurements were analyzed in the Stable Isotope Mass Spectrometry Laboratory at the National Oceanographic Centre, University of Southampton, using a Thermo Fisher Scientific Kiel IV carbonate device coupled to a MAT253 stable isotope ratio mass spectrometer. All samples were measured against the reference standards NBS19 and NBS18, as well as an in-house quality control standard (GS1) and then standardized using a two-point calibration between NBS19 and NBS18 to Vienna Pee Dee Belemnite (VPDB). Long-term analytical precision based on repeat analysis of GS1 is estimated as ± 0.09 for δ 18 O and δ 13 C.

Individual Foraminifera Analysis
To investigate the link between functional traits and environmental change in deep time we use individual foraminifera analysis (IFA) on the extinct, thermocline dwelling, asymbiotic planktonic foraminifera genus Subbotina (Edgar et al., 2013). Subbotina was present globally in Earth's oceans from the early Paleocene through to the end of the Oligocene (∼65-23 Ma) Wade et al., 2011). During the middle Eocene, Subbotinids increased in abundance at multiple sites across the globe flourishing over a period detrimental to many other groups (Macleod et al., 1990;Luciani et al., 2010). A possible reason for this survival is a hypothesized adaptable depth ecology, which is suggested in various intervals of the Eocene based on stable isotope measurements (Macleod et al., 1990;Bralower et al., 1995;Wade, 2004;Dutton et al., 2005;Wade and Pearson, 2008;Stap et al., 2010;Arimoto et al., 2020). Their ecology, diversity fluctuations and depth habitat make Subbotina a versatile group to investigate the link between functional traits and the environment.

Individual Morphological Analysis
To collect morphological traits, 50 individuals of Subbotina were picked from each of the six sample residues (discarding specimens showing gametogenic overgrowth and evidence of reworking) resulting in 300 individuals in total. To obtain 50 Subbotina individuals per sample, each sample was split using a microsplitter until approximately 300 foraminifera remained, shown to represent the diversity of a sample (Al-Sabouni et al., 2007). Subbotina were then picked from this "split." To avoid biasing because of an uneven distribution of individuals on the picking tray, individuals were picked from square cells on the picking tray chosen by a random number generator until 50 well-preserved Subbotina were picked per sample. Subbotina individuals were mounted on glass slides with the aperture facing upward (umbilical view) using double sided sticky tape in groups of up to 20 for morphological analysis (Brombacher et al., 2017). Images of each block of 20 individuals were taken using a Leica M205C stereo microscope with IC90HD camera illuminated from above and then processed using Image Pro 9.1 Software. Automatic measurements of test area (size) and aspect ratio (shape) (Supplementary Figure 1) were taken using the automated image macro in Image Pro (Supplementary Table 3), which have high reproducibility (Brombacher et al., 2017). To maintain a similar scale across all analyses, size was log transformed and mean-centered around 0.

Single-Specimen Stable Isotope Analysis
In addition to the 18 multi-specimen foraminifera samples analyzed, oxygen and carbon stable isotope ratios were determined for 120 Subbotina individuals. During the morphological analysis, described above, each individual from each sample was assigned a number from one to fifty. For each sample 20 individuals were chosen through computer generated random numbers with no replacement. Samples were cleaned in the same manner as the multi-specimen analysis described in section "Stable isotope analysis". Weights of individuals ranged from 13 to 21 µg with an average weight of 14 µg (Supplementary Table 4). To account for the range of weights, vials were loaded into the Kiel in batches, with each batch consisting of vials of individuals with similar weights so that the settings across the whole batch were optimized for a narrower size fraction. Eight samples out of the analyzed 120 failed to record a measurement due to insufficient weight.

Statistical Methods
Ordinary least squares linear models constructed in the R environment (version 4.0.3; R Core Team, 2020) were applied to investigate the drivers of planktonic foraminifera δ 18 O and δ 13 C values. Separate models were built with δ 18 O and δ 13 C as the dependent variable and a combination of environmental (δ 18 O and Age), ecological (δ 13 C) and morphological (test area and test aspect ratio) traits as the independent variables. Numerous models were constructed with varying degrees of interaction. Models were compared using analysis of variance and the best model chosen based on likelihood ratio tests and Akaike Information Criterion (AIC). Clustering analysis was conducted with a Gaussian finite mixture model using the mclust package (Scrucca et al., 2016) with the best model and thus number of clusters chosen using the Bayesian Information Criterion (BIC). The influence of sample age on morphological and stable isotope variables was tested using one-way analysis of variance (ANOVA).

Paleoceanographic Changes Based on Multi-Specimen Analyses
To analyze our trait data in the context of local MECO paleoceanographic changes, we generated stable isotope records from three genera at each of the sites each with different depth habitats (Figure 2 and Supplementary Table 2). Globigerinatheka show the lowest δ 18 O across the interval compared to Catapsydrax and Subbotina (Figure 2A) (Figure 2A).
The vertical thermal structure of the water column can be assessed by calculating the difference in δ 18 O between surface water dwellers (Globigerinatheka) and deeper dwellers (Subbotina (thermocline) and Catapsydrax (subthermocline): (Figure 3A and Supplementary Table 5). A multiple linear regression using the δ 18 O difference between Globigerinatheka and the two deeper dwelling genera (Subbotina and Catapsydrax) as the response variable showed a significant impact on the habitat differences in oxygen isotope space (p < 0.001) and in the MECO (p < 0.01) interval (Supplementary Table 6). There was a predicted 0.71 ± 0.11 overall decrease in the δ 18 O gradient between Globigerinatheka and the deeper dwelling genera (Subbotina and Catapsydrax) across the study interval compared to the gradient change observed in δ 13 C (Figure 3 and Supplementary Table 6). Additionally, there was a predicted 0.56 ± 0.15 decrease in overall isotopic gradient (both δ 18 O and δ 13 C) during the MECO interval compared to the other time slices, which was primarily due to a reduction δ 18 O gradient between all genera (Figure 3 and Supplementary Table 6). An adjusted R 2 of 0.69 shows that 69% of variation in isotopic differences in Figure 3 can be explained by a model that includes isotope grouping (carbon and oxygen) and interval grouping (MECO or not). The δ 18 O surface−Catapsydrax is between 0.52 and 1.30 for most of the interval except at 40.14 Ma where δ 18 O is 0.01 (Figure 3A).  (Figure 2A).
Globigerinatheka consistently show the highest δ 13 C values of the measured genera throughout the study interval, and their values stay within a narrow δ 13 C range throughout (between 2.15 and 2.92 ; Figure 2B) with the highest δ 13 C value recorded immediately after the MECO. Subbotina and Catapsydrax consistently yield δ 13 C values > 1.00 lower than Globigerinatheka and show their highest δ 13 C values in the MECO and immediately after. Subbotina sits slightly above Catapsydrax in δ 13 C space as expected based on their recorded depth habitat but following the MECO this relationship is reversed ( Figure 2B).

Individual Geochemical Analysis of Subbotina
Individuals of Subbotina show variation around "batch" Subbotina measurements in both δ 18 O and δ 13 C space (Figure 2). However, the batch Subbotina measurements plot within 1 standard deviation of the median of the individual foraminifera analyses (IFA; black circle, Figure 2) in each time slice. This is expected and indicates that IFA analysis draws out intraindividual variation within a genus. The range in IFA δ 18 O values at each time slice, apart from at 38.50 Ma, is consistently >1.00 , with the widest ranges observed prior to the MECO where one to three analyses sit outside of two absolute deviations of the median (Figure 2A). These "extreme" values sit apart from the rest of the measurements in each sample. For example, at 41.31 Ma one positive δ 18 O value leads to the largest range across the interval of 2.44 (Figure 2A). Only at 40.14 and 38.50 Ma do all measurements sit within 2.5 absolute deviations of the median. Although these points are more "extreme" there is no reason that they should be excluded. In all samples some individual Subbotina measurements plot within the same space or above that of the Globigerinatheka batch measurements (Figure 2A). In contrast to δ 18 O, individual measurements of Subbotina in δ 13 C space are overall more spaced out with limited or no clusters. The range of δ 13 C values in each time slice is > 1 throughout with the largest range seen between 40.14-41.31 Ma around the MECO. The exception is at 39.56 Ma where it is 0.70 ( Figure 2B). Unlike δ 18 O, Subbotina individuals are consistently below and separated from batch Globigerinatheka δ 13 C values albeit with a reduction in this separation at the MECO ( Figure 2B).
Cross plots of the individual Subbotina δ 13 C and δ 18 O measurements shows an apparent cluster of points in the middle of the plot (0.50-1.25 in δ 13 C and −0.20-1.20 in δ 18 O space) with higher δ 13 C values primarily from the MECO plotting to the side (Figure 4). To determine whether the MECO data genuinely represent a separate cluster of points, a clustering analysis was conducted using a Gaussian finite mixture model. The analysis revealed a spherical, varying volume model with two clusters fitted the data best with a BIC of -255 (Supplementary Figure 2). This model split the 112 data points into two relatively even clusters (Cluster 1: n = 53, Cluster 2: n = 59, Supplementary Table 7). The cluster classification does not cleanly follow the boundaries of our pre-defined time slices (Figure 4), with all but three data points from the MECO time slice and one from 43.50 Ma within a separate cluster (Cluster 1) along with some individuals from each other time slice. The MECO and oldest time slice cluster together as these two samples are climatically similar, as indicated by δ 18 O values, thus the cluster captures the transient warming interval and the earliest stage of the global cooling trend, respectively.
Although clustering analysis identifies two clusters from the data, there is uncertainty particularly where the clusters are in close contact when δ 18 O is around −1.00 (Figure 4). To investigate this signal further, one-way ANOVAs were fitted on carbon and oxygen separately with age as the groups. There was a statistically detectable difference between time slices in both δ 18 O [F(5,106) = 11.84, p < 0.001] (Supplementary Table 10) and δ 13 C [F(5,106) = 9.84, p < 0.001] (Supplementary Table 12). A post hoc TUKEY HSD test showed this difference to only be significant in δ 18 O between the MECO sample (40.14 Ma) and all other samples except for 43.50 Ma at the beginning of the record (p < 0.001) supporting the clustering analysis (Supplementary Table 11). In δ 13 C, the same post hoc test showed no detectable difference between the start (43.50 Ma) and end (38.50 Ma) of the record, but a supported a difference between the MECO and all pre-MECO samples (p < 0.05) (Supplementary Table 13). This, along with the clustering analysis, strengthens the inference that Subbotina isotopic variation is more nuanced in δ 13 C than δ 18 O, likely through the increased role of biological "vital effects".

Understanding Drivers of Individual Foraminifera Analysis
Morphological traits of test shape (test aspect ratio) and size (mean-centered test area on log scale) were also made on 300 Subbotina individuals including the 120 individuals that were used for geochemical analysis. Considering the full morphological dataset, both size and shape show large ranges in all samples with no clear trend with sample age (Figure 5). A one-way ANOVA showed detectable differences between time slices in both mean shape [F(5,294) = 3.171, p < 0.01] (Supplementary Table 14) and size [F(5,294) = 3.185, p < 0.01] (Supplementary Table 15). Used in isolation these traits give us little information, so we integrated the 120 individuals that had both morphological (size) and stable isotope measurements FIGURE 3 | Stable oxygen (A) and carbon (B) isotope gradients between Globigerinatheka (surface) and Subbotina (thermocline) and Catapsydrax (subthermocline) across the MECO interval. As the line approaches zero the isotopic gradient between the different genera is reduced, with 0 indicating no difference in isotopic signature between the different depth habitats. The vertical blue box indicates the MECO. Table 4) into multiple linear regression models because size has an impact on stable isotope and trace element expression in planktonic foraminifera (Elderfield et al., 2002;Friedrich et al., 2012).

(Supplementary
The best supported δ 13 C model included an interaction between δ 18 O and size with an adjusted R 2 of 0.28 (Supplementary Table 16). Of the predictors included, only age at 40.14 and 39.56 Ma differed detectably from the baseline δ 13 C at 38.50 Ma (p < 0.001) (Supplementary Table 17). At 39.56 and 40.14 Ma (MECO) δ 13 C is predicted to increase by 0.46 ± 0.11 and 0.60 ± 0.13 , respectively, compared to δ 13 C at 38.50 Ma (Supplementary Table 17). Our models predict a 0.10 decrease in Subbotina δ 13 C values per log size unit increase, although this relationship is not significant (Supplementary Table 17). Removal of sample age and oxygen isotopes leaves size as the only explanatory variable results in a positive relationship between δ 13 C and size [0.29 ± 0.15 increase per log(size) increase], but not one that is detectably different from 0 (Supplementary  Table 20). Replacing foraminifer size with weight yielded qualitatively similar results with best fitting model structure remaining the same (Supplementary Tables 21-24). The positive δ 13 C-weight relationship (p < 0.05) was qualitatively consistent with the inferred δ 13 C-size relationship, but weight explained more variation (lower residual sum of squares, Supplementary Tables 17, 22) than size.
Since we measured batch isotopes on foraminifera with three typically distinct depth ecologies (surface, thermocline and subthermocline), we can use the batch δ 18 O values to understand the impact of water-depth temperature changes on Subbotina δ 18 O, a proposed thermocline dweller. We used the same model as above but with the age predictor replaced by surface, thermocline and subthermocline δ 18 O measurements as a proxy for water depth temperature. It is the deeper ocean temperatures (thermocline and subthermocline), rather than surface ocean, that drive the individual foraminifera measurement of Subbotina δ 18 O (p < 0.05) (Supplementary

DISCUSSION
Unpicking ecological signals from paleoceanographic signals is difficult, but worthwhile because of the ability to study transient climatic events in deep time. The most common way to infer functionality in deep time is to look for a correlation between a morphological trait and ecological differences (Eronen et al., 2010). We show complicated morphological and geochemical responses to the Middle Eocene Climatic Optimum (MECO, Figures 2-5). The reduction of surface-deep δ 18 O and δ 13 C (1) hydrological changes causing the changes in "batch" isotopes observed but depth habitat of foraminifera remaining constant through the interval (Figure 6) discussed in Section "Paleoenvironmental Changes"; or (2) ecological changes, with deeper dwelling foraminifera (Catapsydrax and Subbotina) migrating upward in the water column to occupy a similar thermal habitat to that of Globigerinatheka (Figure 6) discussed in Section "Measuring functionality in foraminifera traits." The carbon isotope signature (δ 13 C) of planktonic foraminifera is controlled by biology (Edgar et al., 2017) and there are long established relationships between test size and δ 13 C signal in planktonic foraminifera (Berger et al., 1978;Oppo and Fairbanks, 1989;Spero et al., 1991;Elderfield et al., 2002;Friedrich et al., 2012). The δ 13 C-test size relationship is commonly controlled in stable isotope analysis by using narrow size fractions intended to remove the relationship between δ 13 C and size. Unfortunately, this pre-emptive censoring also removes the correlation needed to infer functionality. In this study, we use a wide size fraction (>180 µm) and detect no relationship between Subbotina δ 13 C and size (∼47,219-173,204 µm 2 ) in either multivariate or univariate models when the presence of explanatory variables that track climatic fluctuations are removed (δ 18 O and sample age). Although only a single case study, the lack of size-δ 13 C relationship here has implications for studying functional traits in deep time where vast climatic changes can occur instantaneously in geological terms. Our discussion focusses on the challenges of interpreting these biologically driven isotopic signatures through periods of substantive environmental change.

Paleoenvironmental Changes
Applying trait based studies on thousand-to-million year time scales requires a thorough understanding of climatic and environmental variables. To provide an environmental dimension to our study, we will focus on paleoenvironmental changes (Figure 6) across the 6 Myr focal interval using the inferred depth ecology of individual planktonic foraminifera alongside the measured δ 13 C and δ 18 O of the genera level FIGURE 5 | Morphological [test size (A) and test aspect ratio (B)] variation of all Subbotina measured in this study (n = 300). Size is on a log scale so has no unit of measure. The box represents the interquartile range whilst the whiskers show 1.5*interquartile range. The black line represents the median. Individual filled symbols represent all individual specimens sampled.
isotopes. Through this comparison, we can infer potential water column thermal and trophic state changes across the MECO, which in concert with more temporally resolved studies can provide important climatic and environmental context for our trait-based analysis.

Water Column Thermal Structure Changes Across the Middle Eocene (δ 18 O)
The Eocene saw a gradual change from hothouse to icehouse climate punctuated by short-lived global warming events, the final globally recognized of which is the MECO (Cappelli et al., 2019). Prior to the MECO, a positive shift in benthic foraminifera δ 18 O values at ∼42 Ma in the Atlantic Ocean and other ocean basins (Coxall et al., 2000;Edgar et al., 2007;Cramwinckel et al., 2018;Westerhold et al., 2020) as well as nannofossil assemblage composition changes toward cool water taxa in the Southern Ocean at ∼41.6 Ma (Villa et al., 2008(Villa et al., , 2013 indicate that global temperatures were lower than in the early Eocene and that the transition from a hothouse to icehouse was already underway. Water column cooling resulted in increasing planktonic foraminifera depth habitat stratification (Figures 2, 3).
In our record from the Newfoundland margin, Catapsydrax and Subbotina show progressive δ 18 O increases of 0.45 and 0.61 , respectively, suggesting up to 2 • C cooling of deep water temperatures in comparison to the surface ocean by 41.44 Ma (Figure 2A). Thermal decoupling of the water column due to global cooling would result in the development of a strong, shallow thermocline and more thermally stratified water column prior to the MECO (43.50 to 41.44 Ma; Figure 2A), as suggested in previous research at this locality (Arimoto et al., 2020). In contrast, sea surface temperature change in the subtropical North Atlantic Ocean appear to be relatively minimal with Globigerinatheka δ 18 O of ∼0.23 across the same interval (Figure 2A) supported by minimal surface ocean changes in the South East Atlantic (Galazzo et al., 2014).  (Figure 2A). This suggests that warming was concentrated at thermocline and subthermocline depths (increasing by ∼7 • C) compared to the surface ocean (∼1 • C increase). This muted surface ocean response is not supported by high resolution stable isotope records at this locality that suggest a 4 • C increase in surface ocean temperature (Arimoto et al., 2020) or most other global studies using various geochemical proxies which suggest between 3-6 • C warming of sea surface temperatures across the MECO (Bohaty and Zachos, 2003;Bohaty et al., 2009;Bijl et al., 2010;Cramwinckel et al., 2018Cramwinckel et al., , 2019. Additional sites in the north Atlantic alongside temperature estimates from a greater range of proxies such as Mg/Ca and clumped isotopes are needed to understand the mismatch between this dataset and that of Arimoto et al. (2020) and any site-specific changes in the thermal structure of the water column.
Peak warming in the MECO was followed by globally rapid cooling (Villa et al., 2008;Bohaty et al., 2009;Arimoto et al., 2020) accompanied by increases in cool water taxa (Luciani et al., 2010;Villa et al., 2013), reduced abundances and eventual extinctions of oligotrophic, shallow mixed layer planktonic foraminifera Acarinina and Morozovelloides (Wade, 2004) and calcareous nannofossils discoasters (Villa et al., 2008). In our record, global cooling is represented by an increase in δ 18 O Globigerinatheka−Catapsydrax of ∼1.00 (Figure 3) re-establishing the pre-MECO planktonic foraminifera depth stratification in δ 18 O seen at 43.50 Ma by 38.50 Ma, indicating abrupt cooling of the thermocline following the MECO. Surprisingly, despite such rapid cooling, Catapsydrax takes ∼2 Myrs to reappear below Subbotina in the water column (Figure 2A), possibly not returning to its position below Subbotina in δ 18 O space until thermocline conditions became more stable at 38.50 Ma (Figure 2A Trophic State Changes Across the Middle Eocene Climatic Optimum (δ 13 C) Inferred trophic state changes are highly site specific (Witkowski et al., 2012(Witkowski et al., , 2014Moebius et al., 2015;Cramwinckel et al., 2019). Our study shows that, in tandem with the increased δ 18 O depth stratification between 43.50 and 41.31 Ma [described in section "Water column thermal structure changes across the middle Eocene (δ18O)"], the δ 13 C stratification also increased over the same pre-MECO period by 0.48 ( δ 13 C Globigerinatheka−Catapsydrax ). This increasing separation between genera is interrupted by a transient reduction in δ 13 C Globigerinatheka−Catapsydrax to 1.00 and δ 13 C Globigerinatheka−Subbotina to 0.81 at 40.15 Ma. These δ 13 C values during the peak of the MECO represent the lowest δ 13 C gradients between the surface and deep-water dwellers recorded in our samples ( Figure 3B). The decreasing gradient is a result of increasing Subbotina and Catapsydrax δ 13 C values, coupled with a slight decline in Globigerinatheka of ∼0.14 ( Figure 2B) suggesting a possible decrease in water column primary productivity at Sites U1408 and U1410.
The interpreted reduction in primary productivity from our record during the peak of the MECO event at 40.15 Ma is supported by another study at the same locality (Arimoto et al., 2020) that observed weakened depth stratification between planktic foraminifera and interpreted a large reduction in planktonic foraminifera accumulation rates as an indicator of primary productivity reduction (Arimoto et al., 2020). In addition, decreased productivity during the MECO has been observed in the Southeast Atlantic (Galazzo et al., 2014). Open ocean mid latitude south Atlantic (Galazzo et al., 2014) and north Atlantic (Arimoto et al., 2020) locations are therefore recording a different trophic signal than observed at other sites across the MECO where primary productivity increases have been proposed based on benthic foraminifera accumulation rates and assemblage structure (Boscolo Galazzo et al., 2013;Moebius et al., 2015), increased deposition of organic rich layers in the Tethys ocean (Luciani et al., 2010;Spofforth et al., 2010), increased diatom flux in the Southern Ocean (Witkowski et al., 2012) and North Atlantic Ocean (Witkowski et al., 2014), as well as shifts in planktonic foraminifera communities toward more eutrophic, opportunists (Luciani et al., 2010). At continental margin sites, these changes are attributed to increased weathering and terrestrial input as a result of global warming across the MECO (Moebius et al., 2015). The open ocean setting of our study sites (Sites U1408 and U1410) mean a similar terrestrial input of nutrients is not plausible, which potentially explains the decreased δ 13 C gradient and inferred productivity decrease we observe. Further studies combining micropalaeontological methods and geochemistry are needed in the North Atlantic to understand whether our results represent a local or regional signal. The described thermal and trophic state fluctuations above show the importance of including environmental change in functional trait-based studies. Whilst our study period and data reflect the globally observed rapid changes of the MECO, we also capture global background cooling (Figure 2). While background climatic changes are not often the focus of paleoecological studies, our low-resolution record shows how major changes to water column structure occur before and after large climatic fluctuations. These changes will have an impact on ecosystems and, as we have shown, on traits in planktonic foraminifera.

Functional Traits in Foraminifera in Deep Time
What Does "Functional" Mean for Foraminifera?
This study demonstrates that measurable morphological and ecological traits can be used to infer responses to abiotic forcing in paleoceanographically dynamic environments. However, inferring the functionality of measurable morphological traits is much harder for extinct than for extant species. Further, assigning functions to traits in deep time is often based on observations of extant taxa and the assumption that the observed functional relationship has not changed. Several studies indicate that this approach may be too simplistic Eronen et al., 2010;Edgar et al., 2013). Presence of algal photosymbionts has been shown to be functional (Bé et al., 1977(Bé et al., , 1981 in some modern planktonic foraminifera species (obligate symbiosis) (Takagi et al., 2019), with the only way of determining obligate symbiosis through direct observation. Other functional traits such as spines are not readily preserved in situ and require SEM images to identify. However, more complex morphological traits that are more tightly related to biogeochemical function (such as pore density) can now be easily measured through technological advances (Bé, 1968;Constandache et al., 2013;Burke et al., 2018). While gross morphology has been hypothesized to control buoyancy (Caromel et al., 2014), mathematical models suggest any potential relationship is weak at best (Caromel et al., 2017) and such features are variable even within a constant laboratory environment (Davis et al., 2020). One clear conclusion from these studies and our results is that simple measures of gross test morphology are not primary controls on organismal function, and that interdisciplinary developments offer promising avenues to extract more biogeochemically relevant signals.

Measuring Functionality in Foraminifera Traits
Implying functionality of foraminifera traits is further complicated by foraminifera trait diversity, analytical protocols and trait plasticity. The genus based approach used in this study likely expands the range of morphological or geochemical values compared to species level analyses. Genera have long been argued to represent biological reality (Mayr, 1942) and analysis at the generic level has advantages and disadvantages (Hendricks et al., 2014). Despite their highly resolved, species-level record, a genera based approach is appropriate for planktonic foraminifera as phenotypic and ecological traits are shared across species and genera resulting in morphogroup and ecogroup classifications .
Despite increasing morphological and geochemical niche breadth by measuring genera rather than species, and measuring almost 8 times as many individuals for morphological as geochemical analysis, we do not detect a morphological response in terms of either test size or shape to either the long-term Eocene cooling trend or the transient MECO. In addition, our analysis found no detectable relationship between size and δ 13 C in Subbotina as expected for an asymbiotic foraminifera. Note though that we obtained statistical significance between specimen weight and δ 13 C (Supplementary Table 22), re-emphasizing the importance of measuring the most relevant trait rather than the easiest to measure. This implies that either size is not a functional trait (assuming a δ 13 C-size correlation infers functionality) or plasticity of Subbotina traits is sufficient to mask any functional relationship.
In contrast, to the lack of δ 13 C-size relationship in adult Subbotina, species hosting dinoflagellate algal photosymbionts do have a positive test size-δ 13 C relationship implying functionality (Berger et al., 1978;Oppo and Fairbanks, 1989;Spero et al., 1991;Elderfield et al., 2002;Friedrich et al., 2012;Edgar et al., 2017). The δ 13 C-size relationship in symbiont bearers is a results of algal preferential uptake of isotopically light carbon. If a correlation implies functionality, size is functional in at least symbiont bearing planktonic foraminifera. To understand this relationship further we propose the need for investigations at the individual level outside of analytical size constraints in symbiotic genera and additional research on whether this relationship extends to asymbiotic genera at the individual level.
Individual analysis will also further our ability to constrain the degree of plasticity in planktonic foraminifera traits and therefore better infer their functionality. The discussions in section "Trophic state changes across the Middle Eocene Climatic Optimum (δ 13 C)" focus on linking δ 13 C changes in planktonic foraminifera to trophic changes in the water column. This one-dimensional view of δ 13 C assumes that the depth ecology of planktonic foraminifera does not change and therefore paleohydrological changes drive stable isotope variations. Yet, studies based on extant foraminifera indicate that depth habitat can vary as a result of season, biogeography and environment (e.g., Jonkers and Kučera, 2017;Kretschmer et al., 2018;Taylor et al., 2018;Chernihovsky et al., 2020) as well as changes in life strategy (Darling et al., 2009). Additionally, studies of foraminifera species and across evolutionary lineages have shown changing depth habitats through evolutionary history (Norris et al., 1993;Coxall et al., 2000;Stewart et al., 2012). The increase in δ 13 C values in Subbotina and Catapsydrax ( δ 13 C: 0.57-0.59 , respectively) across the MECO compared to pre-MECO ( Figure 2B), reducing the overall surface-deep δ 13 C depth gradient, suggests these genera could have migrated up in the water column during the MECO supporting an ecological scenario (Figure 6).
An adaptable depth ecology and generalist life strategy has been suggested to explain batch Subbotina variation in δ 18 O and δ 13 C at several points through the Eocene (Macleod et al., 1990;Bralower et al., 1995;Wade, 2004;Dutton et al., 2005;Wade and Pearson, 2008;Stap et al., 2010;Arimoto et al., 2020). In other studies these changes are most often associated with cooling of surface waters and/or increases in productivity (Macleod et al., 1990). However, with no symbionts and a preference for cooler thermocline waters it is hard to envisage this genus moving to a warmer (= shallower) part of the water column.
Utilizing single specimen analysis, we can explore this conundrum further. Depth habitat hypotheses of Subbotina through the Eocene have previously been based on multispecimen (batch) isotope analyses hiding any interindividual variation (Macleod et al., 1990;Bralower et al., 1995;Wade, 2004;Dutton et al., 2005;Wade and Pearson, 2008;Stap et al., 2010;Arimoto et al., 2020). Assuming that our new individual isotopic measurements reflect a genuine isotopic signature of Subbotina, the wide range of individual isotopic values (Figures 2, 4) with some individuals plotting in the same space as both Globigerinatheka (mixed layer) and Catapsydrax (subthermocline) suggests that Subbotina had a large ecological niche extending from the surface ocean to the subthermocline during the middle Eocene (Figure 2). It is this wide ecological niche that may have aided species' resilience through the MECO, multiple early Paleogene hyperthermal events and long-term Eocene cooling. More stable isotope studies of Subbotina at the individual level are needed to understand how the width of the Subbotina ecological niche changed through time and space during the Eocene and through to their demise in the Oligocene.
In this study, we have demonstrated that planktonic foraminifera are an ideal study organism for trait-based studies and can be integrated with paleoceanographic changes to investigate functional trait changes through climatic perturbations. Using an integrated approach with individual based analyses, we have demonstrated profound changes to ecosystems undergoing a transient global warming event. We did, however, detect no evidence that the measured morphological traits, and their relationship to stable isotopes, imply altered functionality across our time period. Further research is needed to reconcile the true meaning of which traits are functional for planktonic foraminifera in deep time, and how we can detect functional relationships statistically in "go-to" morphological traits such as size. We suggest that to truly investigate functionality in planktonic foraminifera we need to be measuring outside of size fractions, using the whole spectrum of genera diversity that planktonic foraminifera offer and leveraging developments in imaging techniques.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
LK, TE, and SN conceived the study. LK analyzed the data, took the lead on writing the manuscript, and collected the data with help from SB and KE. All authors commented on subsequent drafts and approved submission.

FUNDING
This work was supported by the Natural Environment Research Council awards NE/L002531/1 and NE/P019269/1, as well as the University of Southampton's Institute of Life Sciences (IfLS).

ACKNOWLEDGMENTS
This research used samples provided by the International Ocean Discovery Program (IODP). We would like to thank Bastian Hambach and Megan Wilding at SEAPORT laboratories, University of Southampton located at the National Oceanography Centre in Southampton, United Kingdom for assistance in analyzing samples. We would like to thank Anieke Brombacher for taxonomic guidance during initial sample picking. We would also like to thank Max Holmström for taxonomic and analytical guidance as well as helpful discussions.