University of Birmingham Quantifying the stability of planktic foraminiferal physical niches between the Holocene and Last Glacial Maximum

The application of transfer functions on fossil assemblages to reconstruct past environments is fundamentally based on the assumption of stable environmental niches in both space and time. We quantitatively test this assumption for six dominant planktic foraminiferal species ( Globigerinoides ruber (pink), G. ruber (white), Trilobatus sacculifer , Truncorotalia truncatulinoides , Globigerina bulloides , and Neogloboquadrina pachyderma ) by contrasting reconstructions of species realized and optimum distributions in the modern and during the Last Glacial Maximum (LGM) using an ecological niche model (ENM; MaxEnt) and ordination framework. Global ecological niche models calibrated in the modern ocean have high predictive performance when projected to the LGM for subpolar and polar species, indicating that the environmental niches of these taxa are largely stable at the global scale across this interval. In contrast, ENMs had much poorer predictive performance for the optimal niche of tropical-dwelling species, T. sacculifer and G. ruber (pink). This ﬁ nding is supported by independent metrics of niche margin change, suggesting that niche stability in environmental space was greatest for (sub)polar species, with greatest expansion of the niche observed for tropical species. We ﬁ nd that globally calibrated ENMs showed good predictions of species occurrences globally, whereas models calibrated in either the Paci ﬁ c or Atlantic Oceans only and then projected globally performed less well for T. sacculifer . Our results support the assumption of environmental niche stability over the last ~21,000years for most of our focal planktic foraminiferal species and, thus, the application of transfer function techniques for palaeoenvironmental reconstruction during this interval. However, the lower observed niche stability for (sub)tropical taxa T. sacculifer and G. ruber (pink) suggests that (sub)tropical temperatures could be underestimated in the glacial ocean with the strongest effect in the equatorial Atlantic where both species are found today.


Introduction
The distribution and abundance of organisms can provide valuable quantitative and qualitative proxy data for reconstructing palaeoenvironments . In the marine realm, microfossil assemblages are commonly used to estimate a range of palaeoceanographic and palaeoclimatic variables [Imbrie and Kipp, 1971; Climate: Long-Range Investigation, Mapping, and Prediction (CLIMAP) Project Members, 1976; Kucera et al., 2005a]. Planktic foraminifera are perhaps the most common group employed in this manner and have an excellent fossil record, high abundance in deep-sea sediments, and a global distribution [Hemleben et al., 1989]. Their abundance and geographic ranges are strongly related to surface ocean properties, most notably temperature, but also nutrient availability, water column stratification, and turbidity [Ortiz et al., 1995;Bé and Hamlin, 1967;Bé and Tolderlund, 1971;Hemleben et al., 1989;Morey et al., 2005]. Thus, one means of reconstructing past environments is the application of transfer functions. Here a calibration based on the relationship between organisms and their environment in the modern ocean, using the relative abundance of species contributing to an assemblage, is applied to a fossil assemblage to estimate an environmental variable in the past [Berger, 1969;Imbrie and Kipp, 1971;Ortiz and Mix, 1997]. While predominantly applied to the reconstruction of sea surface temperatures (SSTs) [Imbrie and Kipp, 1971;Kipp, 1976;Kucera et al., 2005a], transfer functions also exist for palaeoproductivity [Ivanova et al., 2003] and thermocline depth [Andreason and Ravelo, 1997]. Successful applications of this approach have focused on the Last Glacial Maximum (LGM;~21,000 years ago) [CLIMAP Project Members, 1976;Telford et al., 2004;Kucera et al., 2005b] but have also been extended to the Pliocene where the uncertainty regarding their validity is much larger [Andersson, 1997;Dowsett et al., 2011;Dowsett et al., 2012]. This is because implicit in the transfer function approach is the assumption that the ecological relationship between a faunal assemblage and the environment remains stable [Kucera and Schönfeld, 2007]. This assumption is supported for the last glacial interglacial WATERSON ET AL.
NICHE STABILITY OF PLANKTIC FORAMINIFERA 74 PUBLICATIONS Paleoceanography RESEARCH ARTICLE

10.1002/2016PA002964
Key Points: • Ecological niche models (ENMs) used to test assumption of planktic foraminifer niche stability 21 kya to present • High environmental niche stability of planktic foraminifera between the Last Glacial Maximum and today • Regional calibrations predict global occurrences less well than global calibrations due to niche underrepresentation or multiple genotypes cycles by similar species' environmental optima identified by comparison of abundance and morphometric data [Schmidt et al., 2003] but has not been rigorously tested in an ecological framework. However, even minor changes in a species niche could generate errors in paleoclimate reconstructions and hence impact resulting model-data comparisons and our assessment of climate change impacts [Schmittner et al., 2011].
In fact, niche changes may help to explain some of the mismatches between different temperature proxies at the same location during the LGM [Kucera et al., 2005a; Multiproxy Approach for the Reconstruction of the Glacial Ocean Surface (MARGO) Project Members, 2009].
Quantification of species ecological niches and estimation of niche differences primarily rely on two approaches, ordination techniques, and ecological niche models (ENMs) [Guisan et al., 2014]. Ordination uses direct observations only and compares the difference in environmental attributes at sites where a species occurs in two different time periods to identify any overlap in environmental space [Broennimann et al., 2012]. ENMs are used to estimate species environmental requirements and can be used to describe the ecological niche by relating species occurrence records to the available environment [Guisan and Zimmermann, 2000]. ENMs have been successfully used to assess the retention of niche-related ecological traits in the past [Martinez-Meyer and Peterson, 2006;Nogués-Bravo et al., 2008;Maguire and Stigall, 2009;Stigall, 2012] and to predict present and future distributions of invasive species Weinmann et al., 2013]. Here we use the ENM algorithm MaxEnt which uses presence-background data and compares the environmental conditions at locations of occurrence records with randomly selected points from a background extent to create maps of environmental suitability [Phillips et al., 2006;Merow et al., 2013]. We use the term "niche" here to refer to the multivariate space of physical environmental variables that best correspond to observed species distributions and the associated distribution of potentially abiotically suitable habitats. Other factors, such as biotic interactions and dispersal ability, which limit the occupation of potential niche space, are not incorporated [Araújo and Peterson, 2012]. We augment this analysis with ordination methods which have been shown to more accurately quantify niche overlap (an indicator of niche stability) than ENMs [Broennimann et al., 2012] while being less able to optimize the contribution of different environmental variables to a species' geographical distribution which we analyze using ENMs [Broennimann et al., 2012;Guisan et al., 2014].
In this study we aim to (1) identify the environmental drivers of foraminifer distributions and (2) quantify niche stability between the modern and the LGM using both ordination and ENM techniques. We focus our study on six species dominating the main biogeographic zones (Trilobatus sacculifer, Globigerinoides ruber (pink), G. ruber (white), Truncorotalia truncatulinoides, Globigerina bulloides, and Neogloboquadrina pachyderma). We use MaxEnt to estimate environmental niches in the modern and project these models to the LGM and use an ordination framework to quantify niche margin change metrics of unfilling, expansion, and stability between the LGM and modern. We further test for the presence of nonanalogous environmental space between the modern and the LGM that could bias our analyses, and assess the sensitivity of ENMs to basin-specific versus global calibrations.

Species Occurrence Data
Species occurrence data for the modern and LGM were taken from the Multiproxy Approach for the Reconstruction of the Glacial Ocean Surface (MARGO) project [Kucera et al., 2005a;Barrows and Juggins, 2005;Chen et al., 2005] as archived on the PANGAEA data repository portal (http://pangaea.de/Projects/ MARGO). Total sampled localities for the modern are Atlantic (n = 1166), Pacific (n = 1468), and Indian Oceans (n = 685) and for the LGM are Atlantic (n = 748), Pacific (n = 265), and Indo-Pacific Oceans (n = 273) ( Figure 1). Samples were not included from the Red Sea and the Mediterranean because during the LGM they were characterized by very different environmental conditions [Thunell and Williams, 1989] that lead to the development of abiotic zones that may bias the data sets [Fenton et al., 2000]. All of the data sets included in MARGO have been carefully prefiltered to exclude dissolution-sensitive data points, and thus, we expect preservation to have little impact on our analyses. Certainly, preliminary ENM analyses support little sensitivity of outputs to preservation as a function of depth-dependent dissolution (Figures S3a and S3b in the supporting information). Occurrence data were treated in two different ways; first to describe the total niche (results in the supplementary information) and second, data were filtered to represent the more informative "optimum niche" (Figure 2). This is important as a shift in a species' optimal niche will impact the reliability of Paleoceanography 10.1002/2016PA002964 transfer function application to fossil assemblages [Imbrie and Kipp, 1971;Fraile et al., 2009;Jonkers and Kučera, 2015]. Raw abundance data were converted to relative abundance and then filtered by using defined abundance thresholds for the focal species to create a subset of data representing the species optima ( Figure S7). Thresholds were selected where there were strong increases in species relative abundance compared to background values. Species' optima were then compared to an independent assessment of the optimal niche based on maximum test body sizes determined by Schmidt et al. [2004Schmidt et al. [ , 2006. Total species occurrences are total modern/ optimal modern: N. pachyderma ( ). Geographically filtered occurrences (one presence per environmental grid cell) were used in all MaxEnt analyses to remove the potential inclusion of duplicate records in time or space. Multiple records from one core increase the presence of the species in the analysis of the realized niche as minor changes in abundance do result in differences between absence or presence. These multiple records are unlikely to have any influence on the optimal niche as the differences are very small. Site water depths range from~50 to 5500 m in each ocean basin [Kucera et al., 2005a;Barrows and Juggins, 2005]. Neogloboquadrina pachyderma historically incorporated two morphotypes based on coiling direction which are now known as separate species [Darling et al., 2006]. Neogloboquadrina pachyderma is the sinistral morphotype, and N. incompta is the dextral type; the former is abundant in cooler polar waters, whereas the latter is adapted to subpolar to temperature environments. Here we focus on N. pachyderma sinistral only.   Species-specific ecological characteristics prior to ENM and ordination analyses ensure the inclusion of physiologically important environmental variables within the model framework. Globigerinoides ruber is a (sub) tropical dinoflagellate symbiont bearing species with a temperature optima of~26°C [Bijma et al., 1990;Schmidt et al., 2003] living in the mixed layer. Genetic analyses indicate a wide genetic diversity within this morphospecies (six distinct genotypes in total), suggesting multiple ecologies [Aurahs et al., 2009[Aurahs et al., , 2011. This is clearly highlighted by distinct differences in the distribution and controls on G. ruber (white) and G. ruber (pink) (Figures 3 and 4) in the modern ocean, which are genetically different species [Aurahs et al., 2011]. Globigerinoides ruber (white) is more cosmopolitan in (sub)tropical waters [Kucera, 2007] with temperature optima in the mixed layer between 18 and 25°C [Hecht, 1976;Lombard et al., 2009] and can also tolerate more oligotrophic waters than G. ruber (pink) [Bé, 1959;Bé and Tolderlund, 1971]. Trilobatus sacculifer is a symbiotic species which occupies the mixed layer and dominates tropical oligotrophic assemblages [Reiss et al., 1980;Bijma and Hemleben, 1994]. It exhibits large morphological variability, e.g., individuals with a terminal sac-like chamber (T. sacculifer) and without sac (T. trilobus); however, genetic analyses indicate the presence of just one genotype in the species [André et al., 2012]. Truncorotalia truncatulinoides is a deep-dwelling asymbiotic species that occupies subtropical, temperate, and subpolar regions [Parker, 1962]. Its subsurface habitat in subtropical and tropical waters versus a closer to mixed layer habitat in colder waters [Mulitza et al., 1997] makes its environmental preferences difficult to ascertain [Lohmann and Schweitzer, 1990]. Truncorotalia truncatulinoides is composed of five genetic species adapted to specific hydrographic conditions that are mainly subdivided on the basis of temperature [Quillévéré et al., 2013]. Globigerina bulloides is living in the mixed layer of colder water masses [Bé and Tolderlund, 1971;Naidu and Malmgren, 1996] and highly productive coastal upwelling zones [Thiede and Jünger, 1992]. Consequently G. bulloides distribution patterns are sensitive to food availability [Ortiz et al., 1995;Schiebel et al., 1997] and exhibit two temperature optima at annual SST of 10.4°C and 26.6°C representing its general temperature habitat but also the dominance in tropical and subtropical upwelling zones [Schmidt et al., 2003]. Neogloboquadrina pachyderma sinistral is a polar species with highest abundances between 0 and 3°C [Kucera, 2007;Lombard et al., 2009] and an optimum temperature of~0.7°C [Schmidt et al., 2003] often residing in the pycnocline if food availability facilitates this [Van Nieuwenhove et al., 2016].

Environmental Data
Environmental data are derived from the UK Met Office Unified Model Hadley Centre Coupled Model version 3 (HadCM3); a fully coupled Atmosphere-Ocean General Circulation Model. The oceanic model component has a 1.25°× 1.25°horizontal resolution with 20 vertical levels. The model was run for preindustrial ("modern") and LGM conditions; an in-depth description is in Gordon et al. [2000] and Cox et al. [2001]. The models were initialized from previous preindustrial and LGM simulations and run for a further 1100 years to spin up surface and intermediate waters. Boundary conditions for LGM simulations (greenhouse gases, orbit, ice sheet, and land-sea mask) are based on the Paleoclimate Modelling Intercomparison Project 3 protocol (http://pmip3. lsce.ipsl.fr/). The LGM ocean was initialized from cold ocean conditions; hence, the deep ocean only experiences minimal temperature drift during model spin up (0.4°C in 1100 years). A bilinear interpolation was applied to convert environmental variables from the GCM to 10 minute resolution for use in ENM analyses. This ensured the use of environmental variables at a resolution that adequately captured the foraminiferal physical niche.
For ENM analyses, environmental variables that best describe the ecological limitations on planktic foraminifera distributions were selected [Bé and Tolderlund, 1971;Hemleben et al., 1989;Rutherford et al., 1999;Morey et al., 2005]. Initially, a large number of variables were considered: downward surface solar flux (W m À2 ), mixed layer depth (between 0 and 300 m water depth), potential temperature difference (between 0 and 200 m water depth), brunt-vaisala frequency (measure of temperature stratification in the first 96 m), annual mean temperature (over surface 5 m), temperature seasonality (maximum-minimum monthly temperature), mean temperatures of the wettest and driest quarters, and mean temperatures of the warmest and coldest quarters (over surface 5 m). However, intercorrelated variables are problematic for ENMs, leading to over-fit models and a reduction in predictive ability . Thus, we only retained variables with a Pearson's pairwise correlation coefficient <0.7 (Table S1 in the supporting information). Exploratory MaxEnt analyses found that ENMs incorporating either the mean temperature of the warmest or coldest quarter better predicted LGM foraminifer distributions than annual mean SST. This is unsurprising as absolute Paleoceanography 10.1002/2016PA002964 upper and lower temperature limits more directly impact species and result in lower fecundity, slower growth and smaller sizes [Hecht, 1976;Schmidt et al., 2004;Lombard et al., 2009]. Final ENM analyses therefore use either warm or cold quarter temperatures (in separate calibrations due to high correlations between these two variables) instead of the annual mean, as well as, temperature seasonality, mixed layer depth (MLD), and brunt vaisala frequency. Downward solar flux (W À2 ) was not included because it is highly correlated with SST such that the individual impacts of the two factors cannot be disentangled.
As the carbon cycle is not sufficiently resolved in HadCM3, we were unable to obtain a measure of oceanic net primary productivity for use in the ENM despite foraminiferal abundance and distributions being influenced by food availability [Rutherford et al., 1999;Bé and Hamlin, 1967;Bé and Tolderlund, 1971;Hemleben et al., 1989;Morey et al., 2005]. We considered using the UVic Earth System Climate Model [Weaver et al., 2001;Meissner et al., 2003] rather than HadCM3 in our analyses because it generates both SST and carbonate chemistry in the oceanic component. However, the UVic vertical resolution is too coarse to resolve the surface  Table S8).

MaxEnt
MaxEnt (3.3.3 k [Phillips et al., 2006]) ENMs were calibrated on modern global planktic foraminiferal occurrences and projected to the LGM, using the R "dismo" package [Hijmans et al., 2016]. MaxEnt has shown high predictive performance in comparison with other ENMs, displays good consistency when identifying variable importance, and is a widely used ENM allowing for comparison with other studies [Elith et al., 2006]. A fivefold cross-validation procedure [Elith et al., 2010] was used to create global models for the modern and calculate area under the curve (AUC) statistics (predictive performance measure). AUC values of 1 indicate a perfect model fit, and 0.5 represents a no better than random fit [Fielding and Bell, 1997]. Default MaxEnt model parameters were used, with the exception of the regularization parameter set to 2.0. MaxEnt assesses the importance of each environmental variable to model fit using percentage contribution, permutation importance, and jackknife tests (see Tables S6 and S7 and Figures S1a-S1f) for more details). Geographical extent used for model calibration is typically tied to specific hypotheses of the dispersal capacity of taxa being studied [Anderson and Raza, 2010;Barve et al., 2011]. We applied a global calibration extent to all focal species except G. ruber (pink). This species' modern range is restricted to the Atlantic so training extent was limited to this basin only. A global calibration was chosen for the remaining species on the basis of molecular evidence for genetic mixing of Antarctic and Arctic planktic foraminiferal populations, implying a lack of barriers to dispersal [Darling et al., 2000]. This is supported by fossil studies which indicate that the majority of species ranges are wide-reaching, but large populations are sustained only in regions with suitable abiotic and biotic conditions [Norris, 2000;Sexton and Norris, 1980]. Modern ENMs were projected to global LGM environmental conditions, and a binomial test was used to assess the ability of modern ENMs to predict LGM fossil occurrences. This determines if a model prediction is better than random by assessing the statistical significance of agreement between ENM projections and LGM fossil distributions. Three thresholds were used to define binary environmentally suitable and unsuitable model predictions (0%, or the "least training presence threshold" [in the sense of, Pearson et al., 2007], 90% and 50% fossil occurrences included (see the supporting information for further details).

Testing the Impact of Global Versus Regional Calibrations
The presence of foraminifera genotypes [de Vargas et al., 1999[de Vargas et al., , 2001Kucera and Darling, 2002;Darling and Wade, 2008], which are often ecotypes with a limited distribution compared to the morphotype, may impact transfer function efficiency [Kucera and Darling, 2002]. We tested this by assessing ENM sensitivity to the range of environmental conditions present in individual ocean basins compared with the global range. ENMs were calibrated in the modern Atlantic and Pacific Oceans for Trilobatus sacculifer and Truncorotalia truncatulinoides. These species were selected as extreme cases because T. truncatulinoides has five genetic types [de Vargas et al., 2001[de Vargas et al., , 2004Quillévéré et al., 2013] while T. sacculifer has only one [André et al., 2012]; therefore, any noise introduced by a basin-specific calibration for T. sacculifer cannot be attributed to the presence of cryptic species. Total species occurrences for Pacific calibrations were T. sacculifer (991) and T. truncatulinoides (1037) and for Atlantic calibrations were T. sacculifer (982) and T. truncatulinoides (553). Restricting model calibrations to single ocean basins also impacts the range of the environmental variables incorporated in our ENM analyses. For instance, the Pacific Ocean has a lower salinity and MLD range than the global ocean. Whereas, the Atlantic Ocean is characterized by both lower seasonality and a lower temperature of the warmest quarter than in the global ocean.

Testing Niche Stability: Ordination Framework
Measures of niche margin unfilling and expansion (environmental niche space occurring only in the LGM or modern, respectively) and stability (niche space occurring in both time slices) were quantified by using an ordination framework [Broennimann et al., 2012] with the R "ecospat" package [Broennimann et al., 2016]. These metrics also provide information on the directionality of niche change between the LGM and modern. The framework reduces dimensions in environmental space using a principle components analysis (PCA), and niche quantification analyses are performed within the first two PCA axes. Niche change analyses were calculated relative to the combined modern and LGM species ranges (pooled-range level [Guisan et al., 2014]). Species percentage occupancy per environmental grid cell is calculated by dividing the density of occurrences with the density of the environmental conditions present across the study area [Broennimann et al., 2012]. The framework uses a kernel density function to generate a "smoothed" density of occurrences in environmental space; this helps to minimize unrealistic holes in the niche, which may occur due to low sampling. All ENM and ordination analyses were performed in the free R environment [R Core Team, 2015].
Nonanalogous climates or environments occur over time with climate change, and these can result in unreliable ENM predictions. We tested for the presence of nonanalogue climates using the ExDet software package; this measures the similarity of variables between two time intervals by assessing deviation from the mean and correlation between environmental variables. It also identifies which variables are most important to detected similarity [Mesgaran et al., 2014].

Understanding Foraminifers Modern Environmental Niche
AUC was used as measure of predictive performance for species occurrence data withheld from the model calibration (~20%); average AUCtest scores ranged between 0.78 and 0.90, indicating that the MaxEnt ENMs are able to discriminate presence from background locations   (Figure 3). AUCdiff (difference between AUCtrain and AUCtest) ranged between 0.0006 and 0.0190 (Table S8). For ease we present results for ENMs calibrated with mean temperature of the coldest quarter from herein, unless otherwise stated. Measures of ENM predictive performance (AUC and binomial test scores) and variable importance scores for warm quarter calibrations can be found in Tables S6 and S7 and Figures S1a-f). At the global scale, SST makes the greatest contribution to modern ENM fit for all species (Figure 4). For T. truncatulinoides, while SST was most important (64.7%), seasonality (24.7%) and MLD (10.6%) also make substantial contributions to the model fit.
Models fit at and projected at the regional scale were slightly better or equivalent than species global calibrations, and variable importance was largely similar in the Atlantic and Pacific, for investigated taxa T. truncatulinoides and T. sacculifer. Average AUCtest scores for T. sacculifer were 0.82 for the Atlantic and 0.81 for the Pacific, compared to 0.78 for the global calibration. For T. truncatulinoides, average AUCtest was 0.86 for the Atlantic, 0.89 for the Pacific, and 0.86 for the global calibration. However, regional calibrations showed poorer predictive performance at the global scale for T. sacculifer compared to globally calibrated models (average AUCtest scores for T. sacculifer: 0.52 Atlantic, 0.64 Pacific, and 0.78 global), indicating that the regional calibrations do not capture the global niche as well. In contrast, for T. truncatulinoides the predictive performance of regional calibrations at the global scale was comparable with or better than global calibrations (average AUCtest values were 0.87 (cold) and 0.89 (warm) for the Atlantic, and 0.80 for the Pacific and 0.86 (global) ( Figure 5 and Table S9).

Stability of Foraminifers Niche Through Time
Grid cells with at least one environmental variable outside of the univariate range are confined to small regions close to Antarctica during the LGM ( Figure S10). The brunt vaisala frequency is most influential in these nonanalogue areas ( Figure S5b). Fortunately, very few LGM foraminifera fossil occurrences fall within these nonanalogue regions; therefore, the influence of these areas on the model is expected to be minimal. Modern ENMs were projected to LGM environmental layers to test for stability in the niche between the two time slices. Average binomial test scores were significant at all three thresholds (p < 0.01), for cold and warm temperature quarter calibrations (with the exception of T. truncatulinoides at the least training presence threshold). These results suggest that modern ENMs produce LGM predictions of environmental suitability that are consistent with the distribution of LGM fossil distributions. However, the greatest mismatches between ENM predictions and LGM fossil occurrences are observed for T. truncatulinoides, T. sacculifer, and G. ruber (pink), suggesting that modern models characterize the LGM optimal niche less well for these species (Figure 6).
The cold quarter calibrations show slightly higher thresholds for environmental suitability for tropical species, suggesting that this temperature variable may be more critical in driving niche stability between the LGM and modern for these planktic foraminifera than the temperature of the warmest quarter. Conversely, the warm quarter temperature calibration gives a higher threshold for environmental suitability for N. pachyderma, suggesting that the warmer end of the temperature range better describes the niche of this species in the LGM and modern.
Ordination results suggest that the six planktic foraminiferal species exhibit differing degrees of niche margin stability between the modern and LGM. Neogloboquadrina pachyderma, G. bulloides, G. ruber (white), and G. ruber (pink) have niches more similar than expected by chance (p < 0.05), whereas niche similarity tests were non-significant for T. sacculifer and T. truncatulinoides (i.e., the hypothesis that the LGM and modern niche are no more similar than by chance cannot be rejected) (Figure 7 and Table S10). Species occupying polar and  Table S9).

10.1002/2016PA002964
subpolar waters show the greatest niche stability between the modern and the LGM (N. pachyderma: 99%, G. bulloides: 95%, alongside subtropical species G. ruber (white): 91%) and minimal niche margin expansion between the LGM and modern (0-9%). All species have negligible or no niches occupied solely during the LGM, excluding N. pachyderma, in which niche unfilling is most pronounced (~10%). The largest niche expansion since the LGM is recorded in T. sacculifer (58%), G. ruber (pink) (26%), and T. truncatulinoides (19%). The new habitat is predominantly expressed by a change in niche center in the direction of PC1 (PC1 = 62.35% and PC2 = 11.38%), represented primarily by temperature of the warmest and coldest quarter, and other highly correlated variables including downward solar flux and potential temperature difference ( Figure S6), suggesting that these three species were able to benefit from the increasing SSTs and associated changes in the water column between the LGM and modern (Figure 7 and Table S10). 4. Discussion

Environmental Drivers of Planktic Foraminiferal Distributions in the Modern Ocean
Culture experiments and in situ observations clearly highlight temperature as the most important variable driving modern planktic foraminifer distributions [Bé et al., 1985;Bijma et al., 1990;Ortiz et al., 1995;Morey et al., 2005]. MaxEnt results identify SST as the most important variable on modern planktic foraminifer species distributions (Figure 4), reinforcing the major role this parameter plays in governing biogeographic distributions [Bé and Tolderlund, 1971]. The environmental variables included in our ENM analyses focus solely Figure 7. Planktic foraminifera niche change between the modern and LGM in environmental space determined by ordination techniques. The solid contour lines illustrate the full range (100%) of possible (background) environmental space in the two timeslices, and the dashed lines are 50%. The shading shows the density of species occurrences per grid cell. PC1 accounts for 62.48% total variation and is primarily represented by temperature variables, downward solar flux, and potential temperature difference. PC2 accounts for 11.3% total variation and is represented by the remaining variables: salinity, mixed layer depth, temperature seasonality, and brunt-vaisala frequency (see Figure S6). The blue pixels show the stable niche (common between modern and LGM), the green pixels show the unfilled niche (LGM only), and the red pixels show the expansion of the niche (modern only). The red arrow indicates the change in the direction of the center of the species niche from the LGM to the modern.

10.1002/2016PA002964
on the physical aspects of the foraminiferal niche and omit synergistic factors such as food and light, which modify species distribution under favorable temperatures [Berger, 1969;Bé et al., 1985;Ortiz et al., 1995].
A further potentially important chemical variable that we omit here is seawater carbonate ion concentration, which impacts the geochemical composition of planktic foraminifer tests [Russell et al., 2004;Spero et al., 1997] and calcification [Barker and Elderfield, 2002]. However, its impact on species distribution and diversity is currently unclear. Comprehensive modeling projects are underway which address several of these key oceanic variables [Bopp et al., 2013]; thus, we expect future studies to benefit from those models that incorporate both physical and chemical components of the surface ocean.

Species-Specific Drivers of Biogeography
While temperature is the dominant variable for G. ruber (pink) and G. ruber (white), distinct differences in the distribution of the two genetically different "species" [Aurahs et al., 2011] in the modern ocean are clearly highlighted (Figures 3 and 4). Globigerinoides ruber (white) has a broader temperature optimum (18-25°C [Hecht, 1976;Lombard et al., 2009]) in (sub)tropical waters compared to G. ruber (pink) which has a more tropical distribution [Kucera, 2007]. Our results corroborate the importance of SST as the most dominant variable for T. sacculifer based on experimental results [Bijma et al., 1990]. We find a slightly better model "fit" for T. sacculifer at the basin-specific scale but transferring the calibration to the global scale showed lower predictive performance in comparison to global ENMs, suggesting that geographically-constrained calibration data sets do not capture the entire environmental niche of T. sacculifer ( Figure 5). This is most likely a function of the limited range of environmental space sampled in an individual ocean in comparison to the global ocean. We find that while SST is the dominant variable for T. truncatulinoides, seasonality and MLD also make significant contributions, and is supported by the well-documented dependence of its habitat on water column structure [Lohmann and Schweitzer, 1990;Mulitza et al., 1997;McKenna and Prell, 2004]. Regional calibrations resulted in comparable performance when projected to a global scale for T. truncatulinoides ( Figure 5), suggesting that individual ocean basins encompass the entire environmental range of the morphospecies, as already suggested by the presence of all genotypes in each major basin for this species [Quillévéré et al., 2013]. SST is the dominant variable for G. bulloides model fit, however to a lesser degree than in most other species, as expected due to its dominance in high productivity upwelling areas which explains the relatively high contribution of MLD in our results. Upwelling is highly seasonal at a local scale; thus, this may also explain the contribution of temperature seasonality for G. bulloides. SST is the dominant variable for N. pachyderma model fit, supporting observations that the species' distribution is constrained primarily by its thermal limits [Parker, 1962].
Observational data show that planktic foraminiferal distributions are sensitive to seasonal extremes and that some species respond to environmental change by altering the timing of their presence, with subsequent impacts on predator-prey relationships [Jonkers and Kučera, 2015]. Such responses are important for understanding niche dynamics over time, and we acknowledge that accounting for these is often underestimated by the fallacy of the mean. For those species that are able to alter their depth in response to seasonality changes may still appear to occupy the same temperature despite a shift in the niche.

Ecological Niche Stability Between the LGM and Modern
During the LGM, SST cooling was most pronounced at midlatitudes relative to little or no changes in the tropics and at high latitudes ( Figure S9) [Kucera et al., 2005a].
LGM oceans displayed a higher salinity and greater temperature seasonality in contrast to the modern, with colder winters and an expansion in sea ice extent [de Vernal et al., 2005]. Simulations suggest that the ocean thermal structure changed between the LGM and preindustrial [Telford et al., 2013]. The Western Pacific warm pool contracted, and there was substantive (>3°C) cooling of the eastern equatorial Pacific and eastern boundary currents [Kucera et al., 2005a]. Significant cooling of eastern boundary currents led to steepened SST gradients in the midlatitude North Atlantic, while the subtropical gyres experienced little change.
Our use of optima occurrence data has a substantial impact on ENM predictions to the modern and LGM and ordination analyses, in comparison to occurrence data used in realized niche analyses ( Figures S2-S5). For instance, analysis of species' optimum environmental niche rather than realized niche suggests significant changes in the occupation of physical environmental niche space between the LGM and modern in T. truncatulinoides, T. sacculifer, and G. ruber (pink). This highlights the importance of using species optima for widely Paleoceanography 10.1002/2016PA002964 distributed taxa by removing occurrences that represent marginal niche conditions. Our ENM results suggest that niche stability is best described by cold quarter temperature for tropical species (T. sacculifer and G. ruber) and warm quarter temperature for polar species (N. pachyderma). This has potential consequences for modern analogue techniques as foraminiferal assemblages may represent colder (or warmer) temperatures in the LGM. While this influence is likely to be small it may explain some of the noise observed among proxy signals in some ocean regions during the LGM [MARGO Project Members, 2009].
Based on our findings the optimum environmental niche of three of the six investigated species (N. pachyderma, G. bulloides, and G. ruber (white)) over the past 20 kyr did not change significantly, supporting evidence from biometric [Schmidt et al., 2003] and genetic analyses that the most recent splits in the foraminiferal tree predate the LGM by 100 kyr [de Vargas et al., 2001]. In contrast, ordination results show significant niche expansion for T. truncatulinoides, T. sacculifer, and G. ruber (pink) from the LGM toward the modern primarily driven by temperature. This suggests that these species were able to benefit from increasing SSTs and associated changes in the water column structure between the LGM and present-day, which allowed them to expand their range. ENM predictions to the LGM for T. truncatulinoides ( Figure 6) show fossil occurrences in the equatorial Atlantic in areas of low environmental suitability and might reflect changes in water column structure and thermocline depth [Schneider et al., 1996]. Large areas in the North Atlantic where T. truncatulinoides is commonly found today were unsuitable during the LGM due to changes in the Gulf Stream [Billups et al., 2016].
Modern ENMs for T. sacculifer represent the environmental niche well in tropical regions. Unexpectedly, while the model suggests lower suitability in the South Atlantic subtropical front, an abundance of optima occurrences are found in this area ( Figure S8a). Bé et al. [1981] suggesting that feeding frequency on zooplankton is important for T. sacculifer's optimal niche. Nutrient-rich frontal zones could provide such a habitat; thus, the lack of food availability in our model could help to explain this mismatch.
Modern ENMs are unable to characterize the environmental niche of G. ruber (pink) well or predict its LGM distribution (Figures 3 and 6). This indicates that the variables we use here do not adequately describe the environmental niche of this species. In particular, our analyses indicate high environmental suitability for G. ruber (pink) in the Indian Ocean from which it is absent, indicating that factors other than the physical environment may impact its distribution and that inclusion of biotic factors may help to inform future analyses.
Comparison of GCM reconstructions to multiproxy data sets indicates that the models correctly represent the large-scale features of LGM climate; however, in some areas climate proxies show a wide spread in estimated LGM temperatures. Regions identified by MARGO Project Members [2009] as areas with the largest ranges between different proxies include the equatorial Atlantic and the eastern equatorial upwelling region, which suggest colder temperatures. We observe the greatest niche change in tropical species, T. sacculifer and G. ruber (pink), with expansion from the LGM to the modern, toward regions warmer than those occupied during the LGM. Lower observed niche stability for these taxa means that the absence of these tropical species in LGM assemblages within these regions will likely weight transfer function results toward colder temperature reconstructions.
As you go further back in time differentiation across phylogenies should increasingly display a reduction in conservatism, i.e., niche stability simply breaks down over time [Peterson, 2011]. This will vary between groups and is dependent on the time scales investigated. While our findings suggest that the optimal niche of our focal species is largely (>70%) conservative over the past 20 kyr, the plant record for example, documents nonanalogue floras in the Quaternary [Jackson and Overpeck, 2000;Veloz et al., 2012]. In contrast to the terrestrial realm, there has been limited investigation of niche conservatism in the marine realm with which to compare our results, Saupe et al. [2014] found stability of environmental preferences of marine mollusc species between the Pliocene (~3 Ma) and modern, suggesting that their response to future climate change will be to track suitable habitats or face extinction where the rate of change is too rapid.

Next Steps in Using ENMs in Paleoceanographic Research
A breakdown of niche conservatism is expected over geological timescales longer than the 20 kyr investigated here. Key target intervals for future research are the early Pleistocene, prior to the environmental and evolutionary changes which shaped modern diversity patterns, and past warm interval marine isotope Paleoceanography 10.1002/2016PA002964 stage 11 (424 -374 kyr), which has been the focus of a large data collection and modeling effort; therefore, the necessary environmental and biotic data are readily available for developing ENMs. ENMs have primarily been used with species presence-only data because this is often the most readily available in large-scale biodiversity databases . However, the microfossil community are extremely fortunate in that highly spatially and temporally resolved occurrence and abundance data sets are generated routinely and are available for many taxonomic groups (e.g., coccolithophores, radiolarians, ostracods, diatoms, and foraminifers) unlike in the terrestrial realm where ENMs are most commonly employed. The inclusion of abundance data in ENM analysis would allow the quantification of not just range stability but also how species' optimal niche changes through time (Figure 2), i.e., abundance shifts that may have occurred within the overall niche of a species and could provide vital clues for understanding adaptation and evolution of species.

Conclusions
Global-scale ENM and ordination analyses overall support the assumption of planktic foraminiferal environmental niche stability between the modern and LGM in agreement with previous morphological and genetic approaches. SST is consistently the most important environmental variable to model fit in almost all global ENMs over mixed layer depth, or seasonality this suggests that the thermal limits of the selected species has remained largely constant between the LGM and modern. Ordination results show that the stability of the optimal niche is high for the species that dominate the major subpolar and polar biogeographic zones in the modern ocean. Species occupying tropical regions show lower niche stability and greater niche margin expansion between the LGM and modern. This implies that caution should be taken when applying transfer functions to foraminifera assemblages over this interval that include high abundances of these taxa.
Future studies would benefit from incorporating other factors into our model framework such as light, productivity, and carbonate chemistry. The use of abundance data would also allow us to address range stability and further quantification of the optimal niche through time.