Paleoclimate reconstructions.

The temperature at all three-dimensional locations in prehistoric oceans is not directly observable but can be reconstructed using a coupled AOGCM. We obtained MAT estimates from the UK Met Office Unified Model Hadley Centre Coupled Model (HadCM3), with the MOSES 2.1 land-surface model and top-down representation of interactive foliage and flora including dynamics (TRIFFID) dynamic global vegetation model (52). The precise configuration was HadCM3b-M2.1aD (53). This AOGCM linked atmospheric circulation (at a resolution of 3.75° × 2.5° and 19 levels) to ocean circulation (at 1.25° × 1.25° horizontal resolution and 20 vertical ocean levels). The AOGCM included many environmental variables besides MAT, such as salinity, water age, and N-S and E-W current velocity. Primary productivity and terrestrial nutrient input were not modeled in the AOGCM because these variables are less accurate and more difficult to reconstruct globally, although they are important resource axes for foraminifera in some contexts (36).

A climate simulation was run at 4-ka intervals for the last 700 ka and initialized with climate estimates from similar HadCM3 simulations based on the last glacial–interglacial cycle (54). Each simulation was run for 1,000 y: 900 y of spinup to allow equilibration and a further 100 y that were averaged to produce the output. The ocean layers of interest (within the upper few hundred meters) required only ∼500 y to equilibrate. Deeper waters reached equilibrium more slowly (on the scale of millennia) but impacted surface waters only slightly.

Each simulation at each time step was accorded time-specific boundary conditions for orbital position, greenhouse gas concentrations (carbon dioxide, methane, and nitrous oxide), and ice sheet extent and volume (and corresponding ocean salinity). The greenhouse gas forcing for each time step was linearly interpolated from ice core records. Ice sheet reconstructions were those of De Boer and others (55,56), which are not so well-constrained as more recent reconstructions [e.g., ICE-6G (57) and GLAC-1D (58–60) for the past 27 ka]. However, comparisons of climates at the Last Glacial Maximum (21 ka) using these different ice sheets revealed relatively small changes, except over the ice sheets themselves, and the DeBoer ice sheet reconstructions are the only available estimates of global ice sheets over the last 1 Ma.

Definition of “available environment.”

We defined “available environment” as the set of ocean temperatures from a standard latitude–longitude grid at a given time step, as reconstructed by AOGCM. The first four analyses compared species’ thermal niches to the mean or first difference of mean available environment. Hence, calculation of available conditions affected many sets of results. One issue is that averaging MAT of all AOGCM grid cells indiscriminately in a time step could falsely dampen climatic signals. During glacial intervals, some cells were covered in ice and therefore would not contribute to MAT estimates, bringing up global mean temperature; during warm intervals, some of these polar cells would be open water, bringing down global mean temperature. To obviate this issue of ice-volume fluctuation, we laid a 10° grid of points on the Earth and identified the points that consistently recorded open-water conditions throughout the study interval. We then took the mean of MAT at these ocean points to generate a time series of global temperature, plotted inFig. 3A. For the depth-habitat approach to analyses, we took separate global means at each of the three species’ depth layers.

“Available environment” could be defined as the set of temperatures at all sampled occurrence points instead of at all marine points. Depending on the direction of the null hypothesis, one could argue that using sampled sites would lessen the chance of type I error. In practice, however, we found that replacing the global mean with the sample site mean gave equivalent results. This is because the entire study interval was well sampled (Movie S1), and the mean MAT among all sampled points closely tracked the global mean (Pearson’sr = 0.80), although with slight discrepancies in the timing of some of the interglacial peaks (SI Appendix, Fig. S8). Therefore, results herein are derived from the global marine environment.

Fossil occurrence compilation.

We compiled foraminiferal occurrences from many data sources to span the last 700 ka. Most occurrence records (10⁵) were acquired from datasets in the PANGAEA Open Access library (https://www.pangaea.de/). The compilation also includes records (10³) from the Deep Sea Drilling Project (1968 to 1983), the Ocean Drilling Program (1985 to 2004), the Integrated Ocean Drilling Program (2004 to 2013), and the International Ocean Discovery Program (2013 to present); these were accessed from the Neptune database (https://nsb.mfn-berlin.de/) or the program’s online data repository. Additionally, our dataset includes occurrences from ForCenS, which draws on four compilations and six datasets, including the compilations of climate long-range investigation, mapping, and prediction (CLIMAP), the Brown University Foraminiferal Database, the ATL947 database, and multiproxy approach for the reconstruction of the glacial ocean surface (MARGO) (61).

We harmonized the taxonomy according to currently valid nomenclature. Macro- and microperforate species were included, but benthic foraminifera were excluded. The database contains a text column for “original name” to track synonymization. We excluded two species that originated within the last 700 ka and eight species that went extinct over that span to ensure all study species were extant throughout the entire study interval.

All specimens were collected from sediment cores and core tops, with no live-collection samples, such as from plankton tows. Even the youngest sediment samples are estimated to be centuries or millennia in age (62). Thus, all specimens reflect preindustrial species distributions, which were unaffected by anthropogenic global warming. Paleocoordinates were reconstructed by rotating modern coordinates according to the plate model of Matthews and others (63), in accordance with Neptune methods.

We omitted records from sites at shallow depths where planktonic foraminifera cannot establish viable populations (64). The exact minimum depth for reproduction is unknown and may vary with water conditions. Here, we excluded samples collected at <100 m water depth under the assumption that any foraminifera at those sites drifted in from elsewhere. Deep-water sites could also be problematic because carbonate dissolution could preferentially destroy fossils of certain species. However, such taphonomic artifacts had negligible impact in previous studies of abiotic niche estimates in foraminifera (20), and so we retained deep-water sites for this work.

We extracted MAT values from AOGCM raster cells containing foraminiferal occurrences. If the MAT value of a cell was unavailable (because the cell partially fell on ice or land), we averaged the values of the nine surrounding cells. Finally, we converted presence, absence, and count data to presence-only occurrences for individual species. Absence data were not used subsequently, except to estimate distributions of sampling effort (i.e., sampled core locations as inSI Appendix, Fig. S6 andMovie S1). Occurrences for a species in a time bin were included only if they totaled six or more, a minimum cutoff for robust probability density estimation. Similarly, any species not observed to cross at least six consecutive time bin boundaries was excluded; this cutoff was necessary to ensure sufficient per-species sample size to fit trait evolution models, to correlate time series, and to obtain mean H distances for phylogenetic ANOVA. These sample size constraints reduced the scope of the dataset from 58 to 24 planktonic species. Therefore, results might be more representative of species that are widespread and abundant than those that are restricted and rare.

Spatiotemporal resolution.

Occurrence records were binned into 88 intervals of 8 ka, from the recent subfossil record to 700 ka. The chosen bin length reflects the resolution of the AOGCM outputs (4 ka), the time-averaging of fossil assemblages and imprecision in fossil ages [up to a few thousand years (65)], and the rate of climatic change (fluctuation from glacial minimum to interglacial maximum approximately every 50 ka). To minimize error in bin assignments due to age uncertainty, records were omitted if age estimates were derived from foraminiferal zones or had confidence intervals longer than 2 ka; remaining records could be assigned confidently to a single time bin based on the mean age estimate. Occurrences were paired with AOGCM data from the midpoint of time bins (e.g., climate was modeled at 12 ka for occurrences 16 to 8-ka old).

We rasterized point coordinates to a 1.25° latitude–longitude grid and shortened the dataset to unique combinations of species, grid cells, and time bins. This spatial resolution matched the scale of the AOGCM reconstructions and is appropriate for the magnitude of error in occurrence coordinates. Spatial imprecision could occur from ocean currents displacing foraminiferal tests as the dead organisms sank to the seafloor. For sites at <1 km depth or for large foraminifera, the distance between where the organism died and where it settled may be on the scale of a single grid cell (66,67). For depths of 2 to 3.5 km, a maximum displacement of 100 to 400 km is reasonable, corresponding to an error of up to three grid cells (68,69). Because MAT values are generally quite similar in adjacent cells (i.e., spatially autocorrelated;Fig. 1 andMovie S1), even the largest transport of dead foraminifera to the ocean floor would affect niche estimates undetectably.

Species attributes.

Most living species of foraminifera are benthic, but the ∼50 extant planktonic species are an abundant constituent of the zooplankton at the base of the marine food web. Until recently, all planktonic species were thought to reproduce exclusively sexually, with synchronous release of gametes at the mixed-layer depth of the water column (27). Two research teams have now reported asexual reproduction in planktonic species, but the frequency of asexual reproduction in natural conditions is unknown (70,71). The diversity of ecological strategies in foraminifera includes autotrophy (via association with algal symbionts), herbivory, omnivory, and carnivory (27). In this study, each species was classed into one of five ecotypes (SI Appendix, Table S2), following Aze and others (26): 1) tropical/subtropical, mixed layer, and bearing symbionts; 2) tropical/subtropical, mixed layer, and without symbionts; 3) open ocean and thermocline; 4) open ocean and subthermocline; or 5) high latitude. Ecotype data were used during phylogenetic analysis of variance in niche lability among species.

Foraminifera species inhabit specific and consistent vertical ranges and do not participate in vertical diel migration as do some other plankton, although these protists do migrate to the mixed-layer depth to breed (72). Therefore, we performed each analysis in the following two ways: using MAT at the sea surface to estimate niche differences or using MAT extracted from each focal species’ depth range (a “depth-habitat approach”). In analyses based on the depth-habitat approach, each species’ thermal occupancy was then compared against global ocean temperature from the corresponding depth layer. We searched the literature for quantitative information on each species’ modern depth range to assign taxa to one of the three following depth habitats: surface (40 m in the AOGCM), surface–subsurface (78 m), or subsurface (164 m).SI Appendix, Table S2 lists the depth habitat and reference(s) for each species.

Habitat depth estimates are complicated by the fact that vertical structure of the thermocline, halocline, and photic gradient can shift with latitude, longitude, and time. For instance, a surface-dwelling species might live at 30 m in temperate zones but 10 m in equatorial zones. Our analyses simplified this complexity by using a constant depth point to define each habitat zone. Another assumption of the depth-habitat approach is that species’ depth ranges have been conserved throughout the study interval. Since the goal of the study was to quantify niche lability, there is danger of the following theoretical tautology: assuming stability of depth niches to test the stability of thermal niches. Nonetheless, living depth relates intimately to fundamental aspects of foraminiferal autecology, such as whether a species is constrained to the photic zone by its symbionts, and so depth conservatism on the study timescale is probable. On one hand, the depth-habitat approach incorporates biological complexity into consideration, so results might be more accurate than those based on sea surface data. On the other hand, because of the assumptions of the depth-habitat approach, results from sea surface data might be more precise than those from the depth-habitat approach. Ultimately, the differences between the two methods were small, and temperature values at the sea surface correlated tightly with those from habitat zones (Pearson’sr > 0.95). Previous studies of foraminiferal ecology also recovered equivalent results whether analysis was based on data from the sea surface or at depth (32).

KDE.

We used KDE to construct niche distributions. KDE generates a continuous probability density function (PDF) from a finite set of observations along one or more independent axes. Unlike a probability mass function for discrete data, the PDF at a given point does not represent the probability of an event (e.g., species occurrence) at that point; rather, the integral of the PDF over a defined interval is the probability of an event occurring within those bounds. KDE is a cornerstone of nonparametric statistics with a robust mathematical foundation and applications in many fields, including ecology (37,73–75). KDE can be conducted in any number of dimensions, but every additional dimension requires exponential increases in sample size to maintain low mean square error (73). We chose to prioritize high accuracy of niche estimates (low error) and fine temporal resolution (dividing available data into many time bins) at the expense of a multivariate response (e.g., both MAT and salinity).

One of the difficulties of using KDE with ecological data is the prevalence of nonuniform sampling along environmental gradients. In these instances, one must account for the sampling distribution to make fair and accurate inferences about niches (37,38). This problem is often called “length-bias” in the statistical literature, for the common case in which the probability of observing an event scales linearly with the independent variable (e.g., the chance of identifying an animal increases directly with group size). When using KDE to construct abiotic niches, nonlinear bias could arise from the unequal land cover of different biomes. For instance, because the curvature of the Earth and ocean-atmosphere circulation patterns lead to larger surface areas of temperate than polar habitats, the probability of observing a species in a temperate region of niche space may be higher than that of observing the species in a polar region.

We used the kerneval package in R to construct niches for each species’ occupied MAT in each time bin (version 0.1.2, available athttps://github.com/GwenAntell/kerneval). This tool can correct for uneven sampling using the method of Barmi and Simonoff (76). First, we estimated the PDF of sampling probability along the MAT axis based on the MAT values from all sediment core sites in each time bin. Each core is expected to contain a nearly continuous sedimentary record of ocean conditions and fauna, so we interpolated sampling locations for cores in the middle of their temporal duration for which no records were otherwise present in the database (SI Appendix, Fig. S6A–D). These “pseudo-absences” were used only to estimate the sampling density when standardizing niche estimates; no species presence points were added. Next, the MAT axis for each species was transformed by the inverse of the cumulative distribution function of the sampling probability. We performed KDE for each species in the transformed space and back-transformed the resultant PDF to interpret the niche on the original axis scale.

Choice of bandwidth is a critical step in KDE (73), so we explored densities with rule-of-thumb (nrd0), cross-validation (CV), and Sheather and Jones solve-the-equation (SJ-ste) (77) methods of bandwidth selection. The rule of thumb over-smoothed densities to the point of obscuring differences between species and intervals, while CV generated unrealistically rough PDFs. We choose the SJ-ste method for its middle-ground behavior, although this too could be criticized for undersmoothing densities somewhat (SI Appendix, Fig. S1).

Niche dissimilarity measure.

We used Hellinger’s H distance to quantify the similarity of niche distributions between time bin pairs (Fig. 1 andSI Appendix, Fig. S2). This measure has a long history in the statistics literature, treats the two entities of comparison symmetrically, can be adapted for discrete or multidimensional data, and behaves well even when distributions are highly peaked or contain stretches of zero probability. H ranges from 0 to 1, with 0 indicating identical distributions; this is the opposite interpretation of Schoener’s D, which is a similarity measure rather than a distance measure.

To compare PDFs in a fair way, the sampling universe (observation probability) of each distribution must be the same. In the context of niche overlap calculations, two biological entities must be able to access the same region of niche space for the niche comparisons to be fair (10). If the sampling universe is not shared, one might infer an entity was absent from a region of niche space where another entity was present when, in reality, the region was only available to the latter species. For instance, a species might occupy extremely cold regions during a glacial period but not during an interglacial period merely because there were no equally cold locations on Earth during the interglacial. Without accounting for this difference in availability of environment, one would calculate a larger niche discrepancy than would be justified.

One could standardize niche space either prior to or during niche construction. In either case, one must determine the segment of the niche axis accessible to species in every time bin. The upper and lower limits of the comparison interval are listed inSI Appendix, Fig. S5 for every depth habitat and plotted as dashed lines on niche distributions inFig. 1 andSI Appendix, Fig. S1. If we were to standardize niche space a priori, we would have had to discard any fossil occurrences associated with temperatures outside the interval of comparison (i.e., data points outside the dashed lines inFig. 1 andSI Appendix, Fig. S1). This truncation would remove 19% of presences, almost all of which would exceed the upper limit. Instead, we used all presence points to estimate niche densities but then trimmed the density estimates to a standard MAT interval. We rescaled density functions over this interval, so the probability integrated to unity. Since the standardization took place as the last step of KDE rather than on raw data, time series used in trait evolution models (i.e., the mean and SE of occupied MAT values for each species) were unaffected.

SI Appendix, Fig. S6 outlines the process of calculating the standard temperature interval for niche comparisons. As when estimating the sampling density for niche estimation (KDE), we interpolated gaps in sediment core records with a “range-through” approach (SI Appendix, Fig. S6A–D). We extracted MAT values from AOGCM cells corresponding to core sites in each time bin (SI Appendix, Fig. S6E). We determined the sampled temperature range in each time bin (Movie S1) and then defined the standard interval for niche comparisons as the temperature range accessible to species in all time bins (SI Appendix, Fig. S6F).

For some niche estimation problems, it may be desirable to reflect data across estimation bounds so that the density does not drop toward zero as it approaches those bounds. For example, if the range of sampled values is much narrower than the suspected true breadth of the distribution, then boundary reflection could improve KDE accuracy. We did not reflect boundaries during KDE, however, because the temperature bounds of niche comparison were close to the physiological tolerances of foraminifera. The lower bound for sea-surface data (−1.4 °C) approaches the freezing point of sea ice (−2 °C). The upper bound (27 °C) is slightly lower than the maximum thermal tolerance of foraminifera, 30 to 32 °C (29). However, under default KDE settings, a distribution generally does not fall to zero exactly at each bound, but rather tapers toward zero at a point three times the bandwidth farther away. Thus, for our data, KDE distributions without boundary reflection match closely the laboratory evidence for foraminiferal thermal limits. In addition, 27 °C is the point at which the tight linear relationship between sea surface temperature and thermocline structure begins to break down. Comparing foraminiferal presences to global temperature in a linear way is justified only up to this point (30).

Time series correlation.

The first analysis correlated time series of species’ thermal optima and global temperature. We constructed separate time series for each species, imposing a minimum length requirement of seven time bins. If a species was sampled in multiple, discontinuous sets of at least seven consecutive bins then each segment was analyzed separately to match the methods used in the trait evolution models. For each time step in a species’ series, we used KDE to define the preferred or optimal temperature as the MAT value associated with the PDF maximum. We then correlated each time series of thermal optima with concurrent estimates of mean global MAT. We report the mean and SD of correlation coefficients (Pearson’s r) across all time series (n = 32). A strongly positive correlation coefficient could indicate niche shifts in adjustment to available environment (Methods,Definition of “Available Environment”). In the depth-habitat version of time series analysis, MAT across global marine points was extracted from the same depth layer as was preferred by each species.

Niche lability versus climate change magnitude.

We calculated niche dissimilarity (Hellinger’s H distance) for each species across each boundary it crossed. We then averaged H distances among all species that crossed a given boundary for each of the 87 time bin boundaries. Ultimately, the goal was to relate these niche lability estimates to climate change magnitude, which we defined as the absolute value of the difference in mean global MAT in each pair of consecutive bins. If climate change caused species to alter their niches, one would expect larger niche differences with larger magnitude of environmental change.

We correlated niche dissimilarity and climate change magnitude both with and without removing temporal autocorrelation (“prewhitening”). First, we checked for autocorrelation; an augmented Dickey–Fuller test indicated that the time series of climate change magnitudes was stationary, but the series of mean H distances was not. Therefore, a first-order autoregressive time series model was fitted to the H series; the estimated coefficient of autocorrelation was 0.49. In the case of prewhitening the data to remove autocorrelation, we used the series of residuals from the autoregressive model in place of the original H values as a correlate to climate change magnitude. Temporal autocorrelation was negligible in the time series of H model residuals and the time series of climate change magnitude and in the residuals of a linear model of the former on the latter. This last result confirmed that no further correction was needed to account for pseudoreplication in the correlation of interest.

Niche lability among extreme climate intervals.

We calculated Hellinger’s H distances not only across consecutive time bins but also in pairs of discontinuous time bins representing glacial maxima and minima. We determined the timing of these climatic extremes based on the average global marine temperature curve (Fig. 3A andSI Appendix, Fig. S8 and Table S3). The timings derived from the AOGCM data agreed closely with interglacial onsets reported in the literature (78). For each pair of time bins (e.g., the glacial maximum 28 ka and the interglacial peak 212 ka;Fig. 1), we took the mean H distance among all species observed in both bins. We used boxplots to summarize variation among bin pairs in each of the following three climate categories: glacial–glacial, interglacial–interglacial, and glacial–interglacial intervals (Fig. 3 andSI Appendix, Fig. S9). The mean global temperature difference for time bin pairs in each category was 0.23, 0.14, and 1.99 °C, respectively.

We conducted a variation on the extreme climate comparisons analysis to test sensitivity to sample size deviations through time. One of the benefits of KDE is its invariance to small deviations in sample size; however, the number of observations influences bandwidth calculation to a small degree. The estimate of a complicated distribution will be somewhat rougher with denser sampling than with sparse sampling. This responsiveness to sample size is considered a desirable statistical quality, in theory, because it means that estimates will be perfect representations of distributions as sample size goes to infinity.

Unfortunately, in the empirical case of a comparison of foraminifera niches in the most recent interglacial (4 ka) against any older interval, the difference in sample size was up to two orders of magnitude; of the 42,233 occurrences, 19,711 (47%) were from the 4-ka time step, 2,079 were from 12 ka, and the rest were spread nearly evenly across older time steps. The density estimates for species at 4 ka are noticeably rougher than in much older intervals (SI Appendix, Fig. S1), which could inflate the H distance between densities. In the 86 time bins older than 12 ka, the correlation between sample size and bandwidth estimate was negligible across density estimates (Pearson’sr = 0.040; 95% CI = [−0.0087, 0.089]), and the sample size difference between any consecutive pair of time bins was small enough to be safely ignored in time series and bin-to-bin niche lability analyses. Among glacial intervals, the oldest (668 ka) was sampled least and therefore could influence H distances with other bins. Therefore, we repeated the analysis of niche lability among climate extremes after excluding these two temporal endpoints of glacial–interglacial cycles from comparison (Fig. 3C andSI Appendix, Fig. S9B).

Test for type II error.

We simulated fossil occurrences through time to test for the possibility that the empirical analyses were unable to detect niche changes that occurred (i.e., the possibility of falsely detecting negative signals of niche change, a type II error). In each simulation, species’ occurrence sample sizes and the magnitude of climate change were kept the same as those in the empirical data. We used the geographic occurrences of species at 12 ka as a reference for species’ “true” distributions and assumed geographic ranges were static (and hence, niches adapted to local environmental changes) through time. We chose the time step at 12 ka since global temperatures were close to the long-term mean (Fig. 3A) and sampling was dense throughout species’ geographic distributions.

We generated occurrence data for species in glacial and interglacial intervals (listed inSI Appendix, Table S3), excluding the youngest and oldest extremes (4 and 668 ka) as in the niche lability analysis presented inFig. 3C andSI Appendix, Fig. S9B. For each interval and focal species, we randomly drew the observed number of occurrences for that interval and species, sampling without replacement from the species’ occurrences at 12 ka. Then, we extracted the sea surface MAT at the simulated presences based on AOGCM data for the relevant time step. We used the simulated data to run the “Niche lability among extreme climate intervals” analysis as described inMethods to determine whether mean H distance was larger between pairs of glacial and interglacial intervals than either 1) pairs of glacial intervals or 2) pairs of interglacial intervals. The entire simulation and comparison procedure was repeated 100 times. Boxplot comparisons of results from the first five replicates are illustrated inSI Appendix, Fig. S10A–E.

The simulation framework considers scenarios of complete niche adaptation, but species could respond in many ways beyond a binary of either niche stasis (and geographic habitat tracking) or niche adaptation (and geographic stasis). For instance, species could expand their niche breadth to both occupy the existing geographic range with new environmental conditions and track preferred habitat into new geographic regions. Species could also shrink in realized niche space by occupying only the subset of the original range that remained within tolerance limits. Either of these alternatives would manifest as a difference in the response of the leading and lagging geographic range edge. The simulations conducted here represent a scenario of “lock-step” changes in range edges, which is the most common empirical pattern (46) and allowed us to test whether strong and consistent lability would be detectable with our data and methods.

Trait evolution models.

We fitted five competing models of trait evolution to time series of each species’ MAT occurrences to characterize niche change as static, white noise, wandering, directional, or related to global temperature. We used Akaike weights to compare the five models for each series and defined the “best fit” as the mode of trait change receiving the most support (Fig. 4 andSI Appendix, Fig. S11). Parameters were estimated with a maximum likelihood approach (79,80) implemented in the paleoTS package, version 0.5.2 (81). The functions in this package take the arguments of mean, variance, and sample size of a trait at each time step. These summary statistics must be calculated from the original, finite set of MAT values associated with species’ occurrences; inputting continuous PDFs from KDE into trait models would require substantial rewrites of the paleoTS code. Although optimum temperature (from KDE) arguably has more direct biological meaning than mean temperature among a focal species’ occupied sites, we chose to follow precedent and use mean values with existing code. In practice, the difference between these two metrics is minor; the Pearson correlation coefficient between mean and optimum MAT across all species-bin combinations was 0.85 (95% CI = [0.83, 0.86]) at the sea surface.

The paleoTS models allow several user-specified arguments. We chose not to pool the variance from time steps in a series, since some species’ sample variances were unequal. The singular exception was the covariate-tracking model forGlobigerinoides conglobatus under the depth-habitat approach (marked with an asterisk inSI Appendix, Fig. S11), for which we pooled variances—the specific trait and covariate values for this series by chance resulted in an otherwise unsolvable problem for the default limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm with bound constraints (L-BFGS-B) optimization method. For all series, we specified a restricted maximum likelihood approach (“ancestor–descendant”) instead of a full likelihood approach (“joint”), which treats a series as a joint sample from a multivariate normal distribution. Although the joint approach is more sensitive to directional trends in noisy data (82), package documentation warns that ancestor–descendant parameterization is “generally safer” for covariate-tracking models because it avoids issues of induced autocorrelation (81).

The modeled evolutionary modes were the following: 1) strict stasis, 2) stasis with white noise, 3) unbiased random walk, 4) generalized (i.e., directional) random walk, and 5) tracking of global temperature as a covariate. A strict stasis model assumes that the mean of a trait is constant through time, whereas a general stasis model assumes that the observed mean in each step is drawn from a distribution with a constant mean but positive variance. In other words, in a model of stasissensu latu, observed trait means tend to fluctuate about a constant underlying mean. An unbiased random walk allows the mean in each time bin to move away from the mean in the previous bin by a step length drawn from a distribution with mean of zero and a positive variance; observed trait means migrate from the initial value in a random direction. A general random walk modifies the distribution for the step length such that the mean is constant and nonzero; observed trait means tend to increase or decrease depending on the sign of the step length mean.

For the covariate-tracking model, we identified the series of global mean temperatures during time bins corresponding to the focal species’ series. In the depth-habitat approach, we extracted global temperature values at the preferred depth of the focal species in the respective time bins. The model then estimated a parameter for the correlation between the trait series (species’ occupied temperature) and the covariate series (global temperature). Hence, the covariate-tracking model was akin to the “Time series correlation” analysis, which found species’ preferred temperature to be independent of mean global temperature. Based on this earlier result, the covariate-tracking model might be expected to receive poor support. However, by including covariate tracking in the trait evolution models analysis irrespective of earlier results, we could quantify the relative amount by which other models of niche change were more favored.

It is possible to fit more complicated models than the five considered here. We refrained from fitting punctuations and other complex model types on the grounds of both statistical rigor and relevance to our research questions. Our primary focus was niche change magnitude and covariation with climate, which the five simple trait models already quantify. Incidentally, there is a known bias toward selecting a complicated model when time series are long (83), and the length of time series varied among species in our empirical dataset. If we had included complicated models in comparisons, sampling inconsistencies might have led simple models to be favored more often for rare species and complicated models to be favored more often for common ones.

Phylogenetic ecotype comparisons.

To explore interspecific variation in niche lability, we calculated the mean H across all consecutive bin pairs in which a focal species was observed. We then merged the species-level dataset with the bifurcating morphospecies phylogeny of macroperforate planktonic foraminifera (SI Appendix, Fig. S12) from Aze and others (26), available from the timetree R package, version 3.3.25 (84). Two species were missing from the phylogeny and so were excluded from ecotype comparisons:Globigerinita glutinata andNeogloboquadrina incompta. We performed phylogenetic ANOVA on H values by ecotype (detailed inSpecies Attributes), following the method of Garland and others (85) as implemented in the R package geiger, version 2.0.7 (86). In this approach, the test statistic is calculated in the standard way for ANOVA, while the null distribution is generated by simulating data along the tree. We set the number of simulations at 10,000 to calculate Wilk’s test values. The evolutionary rate was estimated from the mean squared independent contrasts under Brownian motion.