Environmental and biological controls on size-specific 13C and 18O in recent planktonic foraminifera. Paleoceanography

As living organisms, planktonic foraminifera are not passive tracers of the environment. Their test geochemistry — arguably the single most important resource for paleoceanographic research — re ﬂ ects the combined signal of environmental, biological, and preservational processes. For most species, comparisons of test stable isotopic composition within and among taxa provide the primary means for disentangling the relative in ﬂ uences of these different processes. Here we test the foundations of our paleoceanographic interpretations with the ﬁ rst quantitative comparison of the determinants of carbon and oxygen isotopic variation across multiple ocean basins, studies, and species by re-analyzing size-speci ﬁ c data collated from the literature. We ﬁ nd clear evidence of species-speci ﬁ c biological effects (i.e., vital effects), as the intercepts of size-speci ﬁ c carbon and oxygen isotopic compositions differ signi ﬁ cantly among species. Trends in body size and isotopic composition, particularly in dino ﬂ agellate bearing taxa, suggest that much of the size-dependent isotopic variation observed in death assemblages (i.e., core tops and sediments) relates to factors in ﬂ uencing the maximum size obtained by adults rather than ontogeny. The presence and type of photosymbiont hosted (dino ﬂ agellate, chrysophyte, or none) were a major factor affecting species- and size-speci ﬁ c δ 18 O values. In contrast, size-related trends in δ 13 C values were driven by depth habitat (mixed layer, thermocline, subthermocline), symbiont ecology and whether the assemblage was alive or dead when sampled. On this broad geographic and oceanographic scale, ocean basin and biome had a signi ﬁ cant effect on δ 18 O and δ 13 C values . Our analysis and its model-averaged predictions provide a quantitative basis for interpreting size-speci ﬁ c isotopic variation in 22 species of modern macroperforate planktonic foraminifera. We conclude by highlighting existing data gaps and outstanding questions of the relative in ﬂ uence of environmental, preservational, and biological processes on variation in the test geochemistry of planktonic foraminifera.


Introduction
Planktonic foraminifera precipitate calcium carbonate tests. The stable isotopic compositions of these tests are influenced by the ambient environment including the isotopic composition of seawater, temperature, and carbonate ion concentrations [Spero et al., 1997;Spero, 1998]. Because of this, the carbon (δ 13 C) and oxygen (δ 18 O) isotope values of foraminiferal tests are arguably the most important tool available for paleoceanographic reconstructions [Emiliani, 1954;Spero, 1998;Zeebe et al., 2008;Pearson, 2012].
More than half a century of research dedicated to disentangling the determinants of isotopic composition in planktonic foraminifera has led to the great utility of foraminifera in paleoceanography [e.g., Emiliani, Supporting Information: • Readme • Text S1, Figures S1-S28, and Tables S1-S5 • Data set S1 • Data set S2 Correspondence to: T. H. G. Ezard, t.ezard@soton.ac.uk  and Luz, 1982;Jørgensen et al., 1985;Spero et al., 1997] and generally sparse (if any) calibration data for the rest. We also need to account for multiple covarying influences of preservation, oceanography, and biology [e.g., see discussions in Bornemann and Norris, 2007;Birch et al., 2013].
Here we begin addressing this gap by synthesizing existing field data to test whether the empirical evidence supports the purported influence of key hypothesized drivers of isotopic variability in modern planktonic foraminifera. We focus on species-specific size-fractionated isotope curves and interspecific crossplot comparisons (i.e., δ 13 C versus δ 18 O; as in Pearson and Wade [2009], Friedrich et al. [2012], and Birch et al. [2013]) because they are one of the most effective means of constraining the effect of biological and preservational processes on foraminiferal calcite δ 13 C and δ 18 O composition [e.g., D'Hondt and Zachos, 1993;Pearson et al., 1993;Norris, 1996].
To date, modern calibration studies have typically been limited to testing environmental conditions in a single oceanographic location, on a subset of species and parameters [but see Schmidt and Mulitza, 2002].
Here we consider all available data to disentangle the relative influence of environment, biology, and preservation on species-specific trends in size-dependent δ 13 C and δ 18 O relationships. How much evidence is there for key biological characteristics (i.e., depth habitat and symbiont ecology) in size-specific isotopes? How important are environmental and preservational factors in influencing these relationships? We begin by reviewing the background of these existing sets of hypotheses.

Background
Crossplots of carbon and oxygen isotopes from multiple species are used to identify the symbiotic status ( Figure 1a) and relative depth habitats (Figure 1b) of most modern and fossil planktonic foraminifera. This identification relies on the fact that shallow dwelling species typically have relatively low δ 18 O and high δ 13 C values as compared to deeper dwelling species (Figure 1b) Fairbanks, 1992, 1995;Coxall et al., 2007;Aze et al., 2011;Birch et al., 2013]. These interspecific isotopic relationships reflect the isotopic structure of the water column: generally, cooler nutrient-enriched deep waters have higher δ 18 O and lower δ 13 C values relative to warmer surface waters. Vertical gradients in the isotopic composition of seawater generally arise in δ 18 O due to thermal stratification (i.e., fractionation is temperature sensitive: δ 18 O calcite is lower in warm, shallow waters and higher in cool deep waters) [e.g., Pearson, 2012] and in δ 13 C due to the isotopic composition of dissolved inorganic carbon (DIC; δ 13 C of DIC is high in shallow waters where 12 C is bound up in algae and low in deep waters where it has been remineralized) [e.g., Spero et al., 1991].
Correctly identifying foraminiferal depth habitat is paramount for paleoceanographic studies seeking to understand change at a given water depth. However, accurate depth habitat determinations are complicated by biological and preservational factors like species-specific offsets, seasonality, and differential dissolution susceptibility (Figure 1c), which can act to collapse or even invert the order of inferred depth habitat [Savin and Douglas, 1973;Shackleton et al., 1973;Ravelo and Fairbanks, 1995;Pearson et al., 2001;Sexton et al., 2006;Birch et al., 2013;John et al., 2013]. In addition, many biological processes, like metabolic and photosynthetic rates, are directly affected by depth-dependent environmental conditions [Erez and Honjo, 1981;Spero and Williams, 1989;Norris, 1998;Zeebe et al., 2008].
Two well-known biological effects influencing size-specific isotope composition in planktonic foraminiferal species include the formation of gametogenic calcite and the presence or absence of photosymbionts.
Planktonic foraminiferal species with photosynthetic algal symbionts are typically thought to have higher δ 13 C and lower δ 18 O values at any given body size than asymbiotic taxa due to shallower average depth habitats driven by the availability of photosynthetically active radiation (Figure 1a). In addition, species with dinoflagellate endosymbiotic algae typically have steeper size-specific δ 13 C curves, which are attributed to the increasing effect of algal preferential 12 C uptake on the internal inorganic carbon pool with increasing symbiont density or activity through ontogeny [Spero and Deniro, 1987;Norris, 1996Norris, , 1998Bijma et al., 1999;Wolf-Gladrow et al., 1999]. Steep size-δ 13 C relationships are not a foolproof means of identifying photosymbiont bearing taxa: some asymbiotic taxa, including Globigerina bulloides, exhibit relatively large δ 13 C trends (for example, see trends in Spero and Lea [1996] and Birch et al. [2013]); some dinoflagellate-bearing taxa (e.g., Orbulina universa in Birch et al. [2013]) show substantial isotopic variation; and species with chrysophyte endosymbionts are difficult to identify on the basis of their δ 13 C values [Norris, 1998;Bornemann and Norris, 2007].

Stable Isotopic (δ 13 C and δ 18 O) Size Compilation
We compiled measurements of sieve size-specific stable isotope (δ 13 C and δ 18 O) data for 22 extant macroperforate species [e.g., Hemleben et al., 1989;Aze et al., 2011], updating species names in our compilation (Table 1) following the scheme used by Aze et al. [2011]. In all cases, stable isotope measurements were made on multiple individuals. We included all size-specific isotopic data that we could find from the literature, for all modern macroperforate species, as long as three size classes of isotopic data were available. Species-specific justifications are given in supporting information. Note that we distinguish between Globigerinoides ruber pink and G. ruber white (summarized in Darling and Wade [2008]) and the left (sinistral) coiling Neogloboquadrina pachyderma and the right (dextral) coiling Neogloboquadrina incompta [Darling et al., 2000[Darling et al., , 2004[Darling et al., , 2006Darling and Wade, 2008]. Globigerinoides sacculifer (with and without a sac-like final chamber) and Globigerinoides trilobus are all included in G. sacculifer [Hemleben et al., 1987;André et al., 2013].
We examined all modern species with sufficient data (i.e., 22 species), rather than just the handful of species most commonly used in paleoceanography because we aimed to quantify the evidence base underpinning species-specific isotope variability. This necessitated a comparative study across as many species as possible in order to test the effect of shared evolutionary history on observed vital effects with sufficient statistical power. In addition, we included multiple lineages to assess the general applicability of isotopic inferences from modern "workhorse" species to other lineages of foraminifera, as modern taxa are not present in deeper time. Our data compilation is uploaded as supporting information for transparency and to facilitate re-analysis and inspection.
We only considered studies with planktonic foraminiferal isotopes from field-based research that presented three or more sieve-size fractions in order to assess the factors effecting isotopic trends with body size. We deliberately excluded laboratory experiments from the present study as laboratory studies test a portion of the full range of natural environmental conditions and are only available for a small subset of modern foraminifera [Jørgensen et al., 1985;Spero and Lea, 1996;Spero et al., 1997;Bijma et al., 1999;Wolf-Gladrow et al., 1999]. From hereon, we refer to sieve-size fractions as "size fractions" or "size" and recognize that the numerical values do not represent the actual size of a given foraminifera, but rather the midpoint between sieve mesh sizes.

The Database
The database, excluding outliers, features 1358 oxygen and 1147 carbon stable isotope measurements from all the major ocean basins. The data, uploaded as supporting information, include isotopic measurements from recently living individuals (from plankton tows) and long-dead individuals (predominantly from core top samples, including some sediment traps). Fifteen of the 22 species have carbon and oxygen data from both life (plankton nets) and death (traps and cores) assemblages. The influence of this factor (life/death) on δ 13 C and δ 18 O values was tested as a binary explanatory character. The influence of sampling depth, for the plankton tow, sediment trap, or core top, on stable isotopic values was examined as continuous explanatory variables after some preliminary analyses.
Data availability varied substantially among species. For instance, we had from 4 (G. ungulata) to 120 (G. sacculifer) δ 18 O measurements per species (see supporting information and Figure 4 for details). Spatial coverage was similarly skewed: the majority of data come from the Atlantic Ocean, with particularly sparse records from the Caribbean Sea, Pacific Ocean, Southern Ocean, and Weddell Sea (Table 3). Many sample sites are relatively close to land and in regions influenced by seasonal upwelling (Figure 2). The only species sampled in the Caribbean Sea was G. sacculifer and the only species sampled in the Southern Ocean species was N. pachyderma. Since our goal here is a comparison among species and basins, the sparsely sampled Caribbean Sea and Southern Ocean basins were amalgamated with the Atlantic Ocean for the purpose of analyses. "Ocean" basin was therefore considered as a four-way categorical variable contrasting the Atlantic, Pacific, Indian, and Arctic Oceans.
Despite these caveats, the data compilation includes life and death assemblages that are representative of the major lineages of recent planktonic foraminifera, low to high latitude ecosystems, the major ocean basins, and the full range of planktonic foraminiferal ecologies across a range of sampling depths. As such, it provides a means for testing the determinants of size-specific stable isotopic composition in planktonic foraminiferal calcite.
Plotted together, the full data set emphasizes a few fundamental patterns (Figure 3). Small planktonic foraminifera individuals are found throughout the world's oceans, but the largest body sizes only dominate assemblages in the warmer tropical/subtropical biome [Schmidt et al., 2004Al-Sabouni et al., 2007]. The relatively narrow environmental range of the largest individuals is reflected in their narrow δ 18 O range (roughly À2 to 1‰ in the >800 μm fraction) as compared to the smallest body sizes (roughly À3 to 4‰ in the <500 μm fraction) (Figure 3). In contrast to the scattered cone of δ 18 O values, δ 13 C values exhibit a positive, saturating relationship with body size: there is an initial increase in δ 13 C values with body size that eventually plateaus in the largest individuals. We sought to explain the substantial variation in both carbon and oxygen through nonlinear mixed effect models incorporating environmental, preservational, and biological drivers (Table 1).

Model Formulation
We used a statistical modeling approach to derive minimum adequate models (MAMs) and model-averaged predictions of species and size-specific stable carbon and oxygen isotopic trends. Based on existing knowledge of the relationship between isotopic composition and body size (i.e., an increasing linear and a A scarcity of size-dependent oxygen and carbon stable isotope records in the Pacific Ocean, Southern Ocean, and Caribbean Sea is readily apparent. For geographic distribution, see Figure 2. Figure 2. Distribution of available species and size specific stable isotopic measurements for Holocene to Recent planktonic foraminifera used in this study. Core tops and sediment traps (plus) and net tows (cross) plotted against a global bathymetric map extracted from the NOAA server at a resolution of 10 min using the marmap R package (version 0.8). concave-down, nonlinear relationship with size [Berger et al., 1978;e.g., Vincent and Berger, 1981;Birch et al., 2013]) and initial data visualization (Figures 3 and 4), we fitted nonlinear models of the form where y is the isotope of interest (carbon or oxygen) and x is size; a, b, and c are regression parameters that describe the size-dependent change in isotopic composition: a is the intercept, b describes the linear increase, and c the nonlinear (concave down) component of the relationship; ε i is the residual error. Note that we do not force the linear term to increase rather than decrease, nor did we force the nonlinear term to be concave down rather than concave up. Instead, we determined the most likely values of these terms statistically. The subscript i is the species-specific component, which allowed us to test if-and, if so, quantify how-regression parameters varied across the 22 species analyzed. Each biological, ecological, and preservational factor identified (Table 1) was assumed initially to influence each regression parameter (i.e., a, b, and c), with the potential for multiple explanatory factors to influence each parameter. We fitted models using various functions in the nlme [Pinheiro and Bates, 2000] and minpack.lm libraries in the R environment (version 3.0.2 [R Core Team, 2014]). The supporting information contains the data compilation, the annotated computer scripts for model fitting, and the diagnostic plots as justification for our analytical decisions.
All models contain fixed and random effects. Fixed effects are explanatory covariates that influence the mean value of the response variable (i.e., here the environmental, preservational, or biological factors assumed to affect all species); random effects explain the residual variation around those mean trends that cannot be explained by the general fixed effects and therefore is assumed to correspond to (here) differential species-specific responses. Random effects are relevant here because the data are not independent: isotope values within the same species are likely more similar to each other than to isotope values from other species. We aimed to test for species-specific intercepts (parameter a) and linear and nonlinear slopes (parameters b and c) as random effects to determine whether a species (or group of species) differed from the mean trend as modeled by the fixed effects (i.e., environment, biology, and preservation). However, we found that there are insufficient data to adequately test for species-specific size dependence in δ 18 O and δ 13 C values (for full discussion, see supporting information). With the data at hand, we were only able to evaluate species-specific intercepts (i.e., differences in parameter a among species), which were justified from the model diagnostics (see supporting information).
The raw species-specific data ( Figure 4) revealed a great deal of variation in size-specific carbon and oxygen isotope values within and among species. For example, species like Orbulina universa appear to have constrained relationships between size and δ 13 C, whereas Globigerinella siphonifera is more variable. We directly accounted for this difference in variation among species (i.e., heteroscedasticity) and an observed non-independence of model residuals (see supporting information), by extending our basic model formulation to incorporate species-specific variances and autocorrelated errors with size fractions. We could not test for differences among equipment calibrations in different research laboratories.   We fitted models that estimated the influence of both fixed and random effects simultaneously. Initially, the full models for both carbon and for oxygen contained all of the seven biological, ecological, and preservational fixed effects (Table 2) fitted to each of a, b, and c; random species-specific intercepts; species-specific variances; and autocorrelated errors. We then reduced initial model complexity by backward model simplification [Crawley, 2002].
We assess the explanatory power of the candidate environmental, preservational, and biological variables from their contribution to the minimum adequate model (MAM). The MAM was obtained by backward model simplification from the full model using the Akaike Information Criterion (AIC) and likelihood ratio tests [Pinheiro and Bates, 2000;Crawley, 2002]. Here we define the MAM by the lowest AIC score. Testing the global set of models was not possible because numerous combinations do not converge on a solution. The AIC provides a compromise between variance explained and parameters used. Our MAMs for oxygen and carbon (see supporting information) have both the lowest AIC and the fewest parameters.
We found numerous candidate models, which show similar performance in their ability to explain size-specific carbon and oxygen isotope compositions [Burnham and Anderson, 2002, p. 71]. We therefore present model-averaged predictions across all models, weighting by relative likelihood as inferred by the Akaike weights. Akaike weights can be interpreted as the probability that a given model is the "correct" model of those fitted. The model with the lowest AIC has the highest weight. For a more detailed discussion of the statistical models, see the supporting information.

Species-Specific Effects on Size-Specific δ 18 O
Including species-specific intercepts as a random effect improved model fits significantly over models that excluded them, as did autocorrelated residual structures with size and species-specific size-variance relationships (see diagnostic plots in the supporting information). Autocorrelated residual structures with size-specific variance were included because we found that successive size-specific measurements of δ 18 O were more similar than expected by strictly random (white) noise. Similarly, species-specific variance was included by incorporating the different amounts of expressed variation, and therefore residuals, within each of the species in the fitted models (Figure 4). The most influential components of the minimum adequate model (MAM) on changes in δ 18 O values are ocean basin (impacting the intercept and linear slope) followed by sampling depth, live/dead assemblages, biome, and the presence and type of symbionts ( Figure 5). We find no evidence that depth habitat explains a significant amount of the variation in δ 18 O. Two additional models, with different parameter combinations, fell within two AIC values of the minimum adequate model, suggesting they also had "substantial" support [Burnham and Anderson, 2002, p. 71]. The Akaike weight of the minimum adequate model alone is 0.37, indicating just under a 40% probability that this is the "correct" model of those fitted. The three models sum to almost 0.8. The model-averaged predictions  (Table 2), split into contributions to the intercepts (red), linear slopes (blue), and nonlinear (concave) term (green). The bar heights are the difference in Akaike Information Criterion (ΔAIC) between the minimum adequate model (MAM) with and without the focal variable listed: a larger bar indicates that a particular explanatory variable is more influential than a shorter bar. The horizontal grey shaded areas are the criteria of Burnham and Anderson [2002, p. 71] to delineate models: bars within the dark grey region denote a ΔAIC of up to 4 (i.e., model outputs with and without variable are very similar); bars within the grey region denote a ΔAIC of 4-7 ("considerably less" support); the light grey zone is a ΔAIC of 7-10; bars in the white zone beyond the light grey suggest that the model's support is "essentially none," i.e., the focal term explains substantial variation. Variables listed without bars are not included in the MAM. The label for basin in the intercept term is missing because it did not converge on a solution when that term was removed from the MAM.

Species-Specific Effects on Size-Specific δ 13 C
Size-specific carbon isotope trends are more complicated than for oxygen (Figure 7). There is still a crucial role of species-specific random intercepts, differential isotopic variation among species, and autocorrelated errors. As for δ 18 O values, ocean basin (intercept, linear slope) and symbiont type (all parameters) are major determinants of size-specific δ 13 C values. Unlike for δ 18 O values, depth habitat (all parameters) has a substantial influence on the observed variation in size-specific δ 13 C values. Model-averaged coefficients are given in Table 5; the major explanatory driver of size-dependent δ 13 C is the difference between live and dead assemblages (Figures 7 and 8). The summed Akaike weights for the five models with "substantial"  To read this table, note that the intercept includes 0 for continuous effects (water depth) and the first level for categorical variables (i.e., non-spinose, polar/transitional biome, Atlantic Ocean, chrysophyte symbionts, dead assemblages, and mixed layer species). The effect of moving from one of those categories to another is obtained by adding the coefficient, or by adding the coefficient*variable for the continuous variable. Each column added together thereafter gives the parameter value (a, b, c), which you then multiply by size to get predicted size. Therefore, for example, shifting from a chrysophyte-bearing species to an asymbiotic one results in a positive increase in the linear slope parameters, i.e., a steeper linear slope. For species specific intercept adjustments, see Table S2. support is 0.75; the model weight of the MAM is 0.23 indicating less confidence that the model is "correct" than for the analogous δ 18 O analysis. This low confidence in the single "best" model emphasizes the importance of a model-averaging approach that avoids choosing a single model as "correct. In contrast, model-averaged predictions for δ 13 C trends with body size highlight clear differences between life ( Figure 8a) and death (Figures 8b-8d) assemblages. The larger living taxa have a positive, concave-up relationship between body size and δ 13 C. Dinoflagellate bearing taxa have the steepest relationship between body size and δ 13 C, but absolute values generally overlap with other symbiont ecologies at all body sizes (Figure 8a). In death assemblages, δ 13 C values are nonlinear, initially increasing but then plateauing to an approximately zero trend at the largest body sizes (Figures 8b-8d). As with δ 18 O, we find little effect of increasing bottom water depth on δ 13 C trends, as might be expected with increasing diagenetic alteration of foraminiferal tests (Figures 8b-8d).
Model-averaged predictions of reconstructed δ 18 O and δ 13 C values are presented in crossplots in Figure 9 for two scenarios (core top 1500 m and plankton net 250 m) in the Atlantic and Indian Oceans. Two main results are apparent. First, the Indian Ocean predictions yield lower δ 18 O values than their Atlantic counterparts (Figure 9, top versus bottom rows), but the relative order of taxa by symbiont type is consistent.  Table 4. For species specific intercept adjustments, see Table S4.  (Table 2), split into contributions to the intercepts (red), linear slopes (blue), and nonlinear (concave) term (green). The bar heights are the difference in Akaike Information Criterion (ΔAIC) between the minimum adequate model (MAM) with and without the focal variable listed: a larger bar indicates that a particular explanatory variable is more influential than a shorter bar. The horizontal grey shaded areas are the criteria of Burnham and Anderson [2002, p. 71] to delineate models: bars within the dark grey region denote a ΔAIC of up to 4 (very similar); bars within the grey region denote a ΔAIC of 4-7 ("considerably less" support); the light grey zone is a ΔAIC of 7-10; bars in the white zone beyond the light grey suggest that the model's support is "essentially none," i.e., the focal term explains substantial variation. Dinoflagellate-bearing taxa typically have the lowest δ 18 O values in the assemblage, perhaps due to precipitation in warm, well-lit surface waters, followed by chrysophyte bearers, and then asymbiotic taxa with the highest δ 18 O values. Second, the curvature of the species trajectories reflects the nonlinear dependence of δ 13 C and δ 18 O values with size. The deeper prediction (1500 m, core top) has the shallowest within species slopes of δ 13 C to δ 18 O values. Species-specific curves are given in the supporting information.

Discussion
We compiled all available size-specific stable isotope data for modern macroperforate planktonic foraminifera and statistically investigated the controls on test δ 13 C and δ 18 O values across species, space, life history, ontogeny, and death. In summary, our key results are 1. A strong ocean basin effect on both δ 13 C and δ 18 O values (Figures 5, 7, and 9); 2. A clear symbiont effect on size-specific test δ 18 O values (Figures 6 and 9); 3. A significant species-specific effect on δ 13 C and δ 18 O intercepts (parameter a, the basal species offsets; see Tables S2 and S4); 4. An important death assemblage effect modifying apparent δ 13 C to body size relationships and interpretations (Figures 5, 7, and 8); 5. A depth habitat effect on all δ 13 C parameters (Figure 7) but not on δ 18 O parameters ( Figure 5); and 6. The identification of key sampling gaps in the modern calibration set of life and death assemblages ( Figure 2).
The relationships and models presented here are, insofar as is currently possible, an unbiased reflection of the available data across multiple ocean basins on the key drivers of size-dependent changes in planktonic foraminiferal test oxygen and carbon isotopic composition (see the supporting information for fuller discussion and additional figures). While the model-averaged predictions (Tables 4 and 5) can be a used as a predictive tool, we caution against unfettered interpretations because of the inequalities and skew in the available data used to build the models (Figure 3 and Table 3). Our model-averaged results and coefficients should be considered as a set of falsifiable hypotheses, which are rooted in empirical observation and robust statistical methodology, but that are highly dependent upon existing empirical observations. The goal should be to test, reject, and revise these results through targeted field collections and experiments in the future. We structure our discussion around seven emergent questions for ongoing research.

Is Isotopic Variance Driven Primarily by Species-Specific Differences (i.e., Species-Specific Vital Effects) or by General Environmental, Preservational, and Biological Processes?
There is a significant impact of both. All reported models included a species-specific intercept term as a random effect (Tables S2 and S4), as well as several other environmental, preservational, and biological influences (Figures 5 and 7). Species-specific differences in intercepts could reflect underlying biological offsets between taxa [e.g., Erez and Honjo, 1981;Spero and Deniro, 1987;Wolf-Gladrow et al., 1999;Hesse et al., 2014], like different basal metabolic, photosynthetic, and/or calcification rates or other untested environmental and preservational drivers, like seasonality or preservation potential [Berger, 1970;Deuser et al., 1981;Curry et al., 1983;Deuser, 1987;Lohmann, 1995]. It is generally expected that closely related species will be more similar in terms of their biology and ecology, due to their shared evolutionary history, than to more distantly related species [Felsenstein, 1985;Harvey and Pagel, 1991;Freckleton et al., 2002]. Thus, we might expect a significant effect of evolutionary relatedness on the isotopic size-dependency trends in foraminiferal tests that could be used to supplant the "random" species-specific effect. Shared evolutionary history could thus "crack open" the black box of "vital effects" [Spero et al., 1991;Wolf-Gladrow et al., 1999;Zeebe et al., 2008]. Unfortunately, with data only available for 22 species, we had insufficient statistical power to fully test for a role of shared evolutionary history [Pagel, 1999]. The preliminary analyses we could perform suggested a very weak, albeit intriguing, effect of evolutionary history on species-specific intercepts, with much weaker effects on δ 18 O than δ 13 C ( Figure S28 in the supporting information).

Can the Presence of Dinoflagellate and Chrysophyte Photosymbionts be Reliably Inferred From Stable Isotopes?
Yes, with caution. In deep time, the recognition of photosymbiosis relies exclusively on foraminiferal carbon and oxygen isotopes [e.g., Pearson et al., 1993;D'Hondt et al., 1994;Norris, 1998]. Norris [1996] proposed four criteria to identify photosymbiont-bearers in the fossil record: photosymbiont-bearing taxa have (i) the most negative δ 18 O values in the assemblage, (ii) minimal size-related changes in test δ 18 O due to persistent shallower depth habitat throughout ontogeny, (iii) more positive δ 13 C values at all body sizes (again due to shallow depth habitat), and (iv) steeper size-δ 13 C relationships. These criteria are, perhaps, best suited for identifying dinoflagellate bearers, with chrysophyte bearers only weakly (if at all) meeting the criteria [Bornemann and Norris, 2007]. For dinoflagellate-bearing taxa, our results support the cautious use of criteria i Paleoceanography 10.1002/2014PA002735 and iv, as these taxa generally have the most negative δ 18 O values (Figures 7 and 9; particularly true for the largest body sizes) and steepest size-δ 13 C relationships (Figure 8; particularly true for the smallest body sizes). Criteria ii and iii are not supported as useful indicators of photosymbiosis. Asymbiotic and chrysophyte-bearing taxa show less size-dependent δ 18 O change than dinoflagellate bearers ( Figure 6) violating criterion ii. We fail to find consistent, size-independent δ 13 C offsets among the three symbiont groups, violating criterion iii ( Figure 6). For paleoceanographic identification of dinoflagellate bearing taxa, we suggest the combined use of criteria i and iv.
In contrast to the dinoflagellate bearers, chrysophyte-bearing species are hard to distinguish on the basis of their δ 13 C values. Chrysophyte bearers are significantly distinct from asymbiotic species in δ 13 C on all three parameters (equation (1)), but the effect size is so small that it is unlikely to be detected in noisy, empirical, fossil data. Four of the five chrysophyte-hosting species analyzed are thought to be only facultative hosts whereas all five species of dinoflagellate bearing species are obligate photosymbiotic. This difference, compounded by differences in the algal size, density, and location, may affect the carbon subsidy provided by the symbionts to the host during normal chamber formation. In other words, chrysophyte symbionts in planktonic foraminifera may be nearly "invisible" in carbonate δ 13 C because chrysophytes have a lower photosynthetic uptake and thus reduced carbon isotopic fractionation relative to their dinoflagellate counterparts ].
Our analysis suggests, for the first time, that all three groups (dinoflagellate, chrysophyte, and asymbiotic) can be distinguished on the basis of their oxygen isotope values (see characteristic convex δ 18 O curve across ontogeny in Figure 6 and discussion of slopes below). This is surprising, particularly given the relatively weak symbiont influence on δ 13 C values as compared to depth habitat and environmental parameters. Indeed, the presence and type of symbiont is the only biological or ecological factor tested that explains a significant proportion of the variance in oxygen isotope to body size relationship: note the bar heights in Figure 5 and the negative slopes for chrysophyte-bearing species versus neutral (live assemblages) or positive (dead assemblages) slopes for dinoflagellate species in Figure 6. This size-specific oxygen metric for the identification of planktonic foraminifera hosting chrysophyte symbionts is promising but would benefit from further, explicit sampling in the recent and Neogene fossil record.

Can Oceanographic and Ecological Inferences From a Single Site be Applied Globally?
No, particularly not for δ 13 C. Note the impact of basin and biome on δ 18 O, including the missing biome that did not converge (Figures 5 and 7). The stable isotopic composition of seawater is known to vary regionally and by ocean basin Rohling and Cooke, 1999;LeGrande and Schmidt, 2006;Pearson, 2012]. Test isotopic composition is further influenced by local environmental conditions such as the temperature, salinity, and carbonate saturation state of the water from which it precipitates [e.g., Spero et al., 1997;Zeebe, 1999;Hesse et al., 2014]. If the location effect was due solely to seawater composition (e.g., isotopes and carbonate saturation), then it should only impact the intercept term, a, of the size-isotope relationships ( Figures 5 and  7, red bars). For both carbon and oxygen, ocean basin had a significant effect on the linear slope of the relationship (Figures 5 and 7, blue bars), implying that reconstructed size-specific isotopic gradients are, in part, site specific. Multiple mechanisms could underlie site-specific δ 18 O-δ 13 C body size relationships including temperature-dependent metabolic and growth rates [Lombard et al., 2009;Regaudie-De-Gioux and Duarte, 2012;Hesse et al., 2014], carbonate ion effects [Spero et al., 1997;Bijma et al., 1999;Zeebe, 1999], and life history (e.g., seasonality and depth ontogeny) [Sautter and Thunell, 1991;Eguchi et al., 2003;Fallet et al., 2010].
Environment has long been considered to be the primary determinant of δ 18 O and δ 13 C values [e.g., Emiliani, 1954;Hays et al., 1976;Billups and Spero, 1995], as it modifies the external (e.g., seawater composition, temperature, and carbonate saturation) and internal (e.g., respiration and photosynthesis) drivers of carbonate isotopic composition [Wolf-Gladrow et al., 1999;Hesse et al., 2014]. Our results underline the importance of considering environment when using size-specific isotopic relationships in fossil foraminifera, as the slope of this relationship, even in well-preserved samples, can vary.

Does the Difference in δ 18 O and δ 13 C Values Among Species Reflect Their Relative Depth Habitat?
From the current compilation, the evidence is mixed. In our analyses, depth habitat failed to account for significant variation in δ 18 O ( Figure 5) but influenced all parameters for size trends in δ 13 C (Figure 7). The lack of evidence for a depth habitat influence in δ 18 O may arise because the horizontal geographical variation overwhelms the vertical signal, or may reflect an underlying limitation of isotopically defined depth habitats.
Sixty years ago, Cesare Emiliani first inferred differing depth habitat preferences of planktonic foraminifera using the average δ 18 O composition of species and assuming isotopic equilibrium with the surrounding waters [Emiliani, 1954]. Then, as now, stable isotopic measurements provided a convenient means for identifying a single "average" depth habitat for each species, an inference critical for determining depth-specific temperatures and water column stratification [Hemleben et al., 1989;Mulitza et al., 1997Mulitza et al., , 2003. In reality, planktonic foraminiferal populations live across a broad range of water depths that vary between locations, with ontogeny, and with water column structure, among others [Duplessy et al., 1981a;Ravelo and Fairbanks, 1992;Niebler et al., 1999;Cléroux et al., 2013]. Many species have highly dynamic populations, with large differences in population abundance and adult body size across seasons and years [Deuser et al., 1981;Eguchi et al., 2003;Mohtadi et al., 2009]. In our study, as in all strictly paleoceanographic calibrations, seasonal signals are unaccounted for but could isotopically overwhelm the effect of depth habitat. In short, our failure to find a significant depth habitat may arise because oxygen isotopes are less informative of depth habitat than has long been hypothesized (a criticism first raised by Shackleton et al. [1973]). The alternatives are (1) that our data are simply insufficient to test the effectiveness of δ 18 O values for assessing foraminiferal relative vertical position in the water column or (2) that by subdividing depth habitat into three general categories, mixed layer, thermocline, and subthermocline, we obfuscated fine, continuous, vertical differences in depth habitats recorded in δ 18 O. More data are required to assess the role of depth habitat on species' size-specific δ 18 O through paired tow, sediment trap, and core top studies.

Is There a Diagenetic Effect on Test δ 13 C and δ 18 O Values?
Probably not in the core top samples we examined. Our model-averaged predictions indicate that the overall pattern of test size-δ 13 C and δ 18 O values and the depth ordering of species are consistent across the range of tested water depths (1500-4500 m; Figures 6 and 8). In particular, δ 18 O-size gradients in dinoflagellate-bearing species are steeper (particularly at small sample sizes) than those of chrysophyte-bearing species, regardless of water depth ( Figure 6). However, our results and inferences are notably restricted to relatively well preserved core top sediments and not the more extensively altered specimens typically found in sediment cores. Furthermore, our focus is on a comparison among ocean basins and between living and recently dead assemblages, rather than the effects of dissolution below the lysocline or in sediments.
As carbonate preservation varies strongly with water depth and geography due to variations in carbonate saturation state and sedimentary conditions [e.g., Archer, 1996;Ridgwell, 2007], and both have a significant effect on δ 18 O trends ( Figures 5, 6, and 9), we briefly consider how depth-dependent isotopic changes might be attributed to diagenetic factors. Postmortem dissolution of planktonic foraminiferal tests in the water column and at the seafloor should lead to progressively higher δ 18 O and lower δ 13 C absolute values with increasing water depth [Savin and Douglas, 1973;Berger and Killingley, 1977;Bonneau et al., 1980;Wu and Berger, 1989], with smaller and/or thinner-walled taxa being more susceptible than larger ones [Berger and Piper, 1972;Berger et al., 1982;Wu and Berger, 1989;Hönisch and Hemming, 2004]. In contrast, replacement of planktonic foraminiferal calcite by inorganic calcite precipitated in situ at the seafloor (recrystallization) will result in overall higher δ 18 O and δ 13 C values with the most pronounced impact on δ 18 O values [e.g., Schrag et al., 1995;Pearson et al., 2001]. We might also expect that the surface-dwelling, symbiont-bearing taxa show a more pronounced change with increased bottom water depth than the deeper-dwelling, thicker-walled taxa such as G. tumida or P. obliquiloculata. Although we do observe a significant effect of water depth on δ 18 O intercepts (parameter a; Figure 5), the effect is primarily driven by differences in the water depth of living, recently dead, and long-dead assemblages (i.e., there is a noticeable difference between Figures 6a and 6b-6d but not among Figures 6b-6d). Thus, we fail to readily observe expected diagenetic effects and, as such, we cannot easily explain the patterns observed in the dataset by differences in the degree of dissolution or recrystallization at the seafloor.

Does the Different Composition of Life and Death Assemblages Change the Underlying Drivers or Foraminiferal Stable Isotopic Trends?
Most likely. Here we find that the contrast between living and dead assemblages is the most important factor influencing the relationship between body size and δ 13 C (i.e., parameters b and c; Figure 7) and one of two factors affecting basal δ 18 O values (parameter a; Figure 5). We suspect that the importance of life versus death assemblages arises because ontogenetic factors like growth rate, symbiont density, and the addition of gametogenic calcite structure size-isotope relationships in living assemblages. In contrast, other factors structure size-isotope relationships in death assemblages, namely those that effect the ultimate adult body size attained by individuals. We dub this difference in isotopic composition, and the factors influencing them, between life and death assemblages the "Spero life-death effect" (SLDE).
In living assemblages of foraminifera sampled by nets or tows, the vast majority of individuals are small (<125 μm) [Emiliani, 1971;Ravelo and Fairbanks, 1995], and variation in size correlates to ontogenetic stage. In other words, small individuals are small because they are young and large individuals are large because they are old. The opposite is true in sedimentary death assemblages, where the majority of planktonic foraminifera are large (>125 μm) [Berger, 1971;Bé et al., 1985]. In these fossil assemblages, much of the observed size variation occurs in postgametogenic individuals. Medium to smaller individuals are small because they were small adults at the time they reproduced. Large individuals are large because they were large adults when they reproduced [e.g., Berger, 1971].
Despite known differences in the composition of life and death assemblages, size-related stable isotopic trends in fossils are often extrapolated to living populations using ontogenetic hypotheses such as developmental shifts in the relative importance of metabolic effects versus photosymbionts, ontogenetic depth habitat variation, and the addition of gametogenic calcite [e.g., Schweitzer and Lohmann, 1991;Norris, 1998;Birch et al., 2013]. The importance of life versus death assemblages to predicting size-specific isotopic composition (Figures 5 and 7) supports the consideration of other mechanisms for size-specific trends in death assemblages.
The size-related δ 18 O and δ 13 C trends in dinoflagellate bearing species (Figures 6 and 7) follow the pattern predicted by Spero and others [Spero and Williams, 1988;Spero and Lea, 1993] for isotopic variation among same-stage adults. By varying the amount of incident photosynthetically active radiation in lab cultures, Spero and Williams [1988] found they could alter both the maximum test size attained and the extent of 13 C enrichment in the dinoflagellate-bearing species Orbulina universa. From this, they reasoned that, in death assemblages in nature, the largest adults should have high δ 13 C and low δ 18 O values due to growing in high irradiance and relatively warm surface oceans as compared to the smaller, individuals with low δ 13 C and high δ 18 O values growing closer to the thermocline [Spero and Williams, 1988]. Other work has also shown that variation in adult body size is driven by conditions experienced over the lifespan of the individuals [Caron et al., 1987a[Caron et al., , 1987bLombard et al., 2009], that adult size relates to symbiont presence and activity [Caron et al., 1981;Bé et al., 1982;Spero and Lea, 1993], and that similar relationships between irradiance, body size, and isotopic composition occur in other taxa (e.g., Globigerinoides sacculifer and Globigerina bulloides in Spero [1992], Spero and Lea [1993], and Bemis et al. [1998]). These mechanisms, as detailed by Spero and others, are not ontogenetic mechanisms and the patterns observed cannot be attributed to the addition of gametogenic calcite alone. Rather, we consider the factors explored by Spero and others among those that affect the size and isotopic composition of same stage adults-that is, they are the mechanisms behind Spero life-death effects.
Spero life-death effects have already been shown to account for size-specific variation in δ 13 C in downcore assemblages (e.g., G. sacculifer as demonstrated in Spero et al. [2003]), with critical implications for paleoceanographic reconstructions. In our data, the Spero life-death effect is particularly evident in dinoflagellate-bearing taxa, which exhibit a positive δ 13 C to body size ( Figure 8) and negative δ 18 O to body size ( Figure 6) relationship. For example, compare the saturation in δ 13 C with increasing body size in death assemblages (Figures 8b-8d) to the apparently unbounded increase in living assemblages (Figure 8a).
The Spero effect directly contrasts with the ontogenetic hypothesis for trends in stable isotope composition with body size, which predicts-among other things-that the largest adults should have relatively high δ 18 O values due to reproduction (and the addition of gametogenic calcite) at depth. Ontogenetic hypotheses implicitly attribute the variation in body size, even in death assemblages, to differences in life history stage, i.e., the age of the individual foraminifera. Were the patterns in our data principally explained by ontogeny, we should be able to observe two distinct body size-δ 13 C, and δ 18 O relationships: (i) a positive relationship at small body sizes due to metabolic and kinetic effects and (ii) a positive relationship between each of δ 13 C and δ 18 O and the largest body sizes due to increasing symbiont density/activity and deeper habitat Paleoceanography 10.1002/2014PA002735 preference in the water column. Regardless of symbiont type, we do not observe either of these relationships in oxygen ( Figure 6) and only the former in carbon (Figure 8).
The importance of Spero life-death effects does not supplant all ontogenetic mechanisms, rather it emphasizes an additional set of factors to consider, particularly in individuals greater than roughly 200 μm. For instance, Schweitzer and Lohmann [1991] found a mix of pregametogenic and postgametogenic individuals when examining M. menardii and G. tumida less than~300 μm. In addition, juvenile foraminifera are common in the smaller size classes (<125 μm), but these are rarely measured for isotopic composition due to the difficulty of accurate identifying species and obtaining sufficient material for analysis (as discussed in Friedrich et al. [2012]).
The Spero life-death effect highlights an important issue for investigations of seasonality and interannual variation. In most paleoceanographic studies [e.g., Killingley et al., 1981;Koutavas et al., 2006;Ganssen et al., 2011;Khider et al., 2011;Sadekov et al., 2013], long-term dynamics in seasonality and/or interannual variation is considered by examining isotopic variation through time within a narrow size, or weight, class of dead individuals. However, if variation in adult body size is primarily driven by environmental conditions, such studies exclude the very individuals that record the dynamics the authors seek to infer.
There are several means for identifying the relative importance and drivers of Spero life-death effects and ontogenetic mechanisms in future studies. These include comparing the size-specific isotopic trends in depth-stratified plankton tows and sediment traps, the use of single chamber laser ablation techniques (as in Houston et al. [1999]), the separation and independent analysis of postgametogenic (for SLDE) and pregametogenic (for ontogeny) individuals using shell weights (as in Schweitzer and Lohmann [1991]), and the careful consideration of temporal averaging across seasons, years and centuries (for examples of refined seasonal approaches, see Fallet et al. [2010]).

Are There Clear Targets and Gaps in the Existing Data Set That Would Allow Us to Refine Stable Isotopic Interpretations?
Yes. Although all the major ocean basins are represented, sampling is highly biased to death assemblages from shelf and upwelling environments in the Atlantic and Indian Oceans (Figure 2 and Table 3), leaving open ocean gyres alarmingly undersampled. In addition, there are relatively few studies that directly compare size-specific isotopic trends in tow data with sediment traps and core tops (important direct comparisons include Williams et al. [1981], Bouvier-Soumagnac and Duplessy [1985], and Ravelo and Fairbanks [1992]). With limited paired life and death assemblage data, the importance of depth habitat occupancy (i.e., vertical niche partitioning) and preservational factors as they influence size-specific isotopic trends is largely untested. It is unclear from the existing net tow literature whether species subdivide the water column as often inferred from stable isotopic data. Similarly, a more direct study of environmental effects including temperature, carbonate saturation, salinity, and nutrients is needed. With the current data compilation, we had the statistical power to test for coarse environmental effects through biome, basin, water depth, and depth habitat (Table 1).
For many species, there is a paucity of data (e.g., G. rubescens, G. ungulata, and P. obliquiloculata; Figure 4) or a high degree of scatter (e.g., G. bulloides, N. incompta, and N. pachyderma; Figure 4), further limiting our ability to make general inferences. The available size-specific isotopic data are predominantly from core top sediments (supporting information), but our understanding of the various controls on isotopic fractionation in calcite is developed primarily through ontogenetic studies of live individuals in culture and target species in nets. This biases our interpretation of death assemblages toward ontogenetic hypotheses rather than Spero life-death effects. In order to document and evaluate both ontogenetic and Spero life-death effects, size-specific calibrations across an expanded size range (as per Friedrich et al. [2012] and Birch et al. [2013]) of paired living and death assemblages from open ocean gyres and underrepresented ocean basins are needed.
In contrast, understanding the evolution of species-specific vital effects requires additional sampling of unrepresented species in the Atlantic, to provide a sufficiently dense record in a single basin (Table 3).
Our supporting information includes the data compilation and annotated computer code to recreate the models fitted and discussed here. We hope that this resource encourages future studies to integrate new data into the wider context provided here, so that our predictions can be re-evaluated as additional data become available.

Conclusions
Modern planktonic foraminifera provide a ready means for generating, testing, and benchmarking theories used to infer paleoclimatic and paleoecologic change. Past studies have typically focused on the generation of new data for a few species or regions. Here we compiled all available existing data, synthesized existing knowledge using integrated statistical modeling, and examined the evidence base that underpins our understanding of isotope geochemistry in this important study system. Our results suggest that size-dependent δ 13 C and δ 18 O trends are mediated predominantly by differences among oceanic regions (basin and biome), the presence and type of symbionts, and for δ 13 C, depth habitat, sample water depth and whether the assemblage was alive or dead. We find compelling evidence, dubbed here the Spero life-death effect, for the dominance and importance of death assemblage factors in shaping size-dependent stable isotopic curves in fossil assemblages. Species-specific factors have a significant effect on the basal δ 13 C and δ 18 O values of taxa. Rather surprisingly, symbiont ecology is more important than depth habitat in explaining δ 18 O values. δ 18 O-size trends may represent a new tool for identifying the presence of chrysophyte-bearing taxa in the fossil record. Depth habitat is at least as important on δ 13 C values than symbiont ecology, which also contradicts typical expectations.