The Dynamics of Diachronous Extinction Associated With Climatic Deterioration Near the Neogene/Quaternary Boundary

To predict extinction, we must understand the processes leading to terminal population decline. Once a critical threshold of population size is reached, small environmental perturbations can push a species over the cliff‐edge to extinction, so the main drivers of extinction are the factors that cause the initial reduction in population size. Most studies of population decline leading up to extinction focus on modern species in a human‐dominated world. The drivers of population decline leading to non‐human mediated extinctions are less well known but changes in climate are arguably the most widely invoked mechanism. Here, we report data on >16,000 individuals of the planktonic foraminifer Globoconella puncticulata from six sites in the Atlantic Ocean along an 83 degree‐long latitudinal transect, over a 600,000‐years interval leading up to the species’ global extinction during the late Pliocene‐earliest Pleistocene intensification of Northern Hemisphere glaciation. We show changes in geographic range, abundance, and body size. We find that populations do not follow a North‐to‐South sequence in extinction as Earth cooled and developed large ice sheets in the high latitudes of the Northern Hemisphere. Instead, our results suggest that (a) populations are differentially adapted to local environmental conditions such as nutrient availability, (b) population dynamics in core populations differ from those at the edge of their range, and (c) individual population responses to external pressures are essential to understanding the drivers of global extinction. Our study demonstrates the potential to transform our understanding of extinction dynamics through spatially replicated sampling of the highly resolved marine microfossil record.

the point where small external perturbations can push a species over the cliff-edge to extinction (Gilpin & Soulé, 1986;Hedrick et al., 1996;Raup, 1981). Because these modest perturbations would likely have been more readily absorbed by healthy populations, the main drivers of extinction are the factors that reduce population size initially and thus enhance extinction risk (Hedrick et al., 1996;Raup, 1981). Therefore, to understand extinction, long-term records of biotic and abiotic change preceding extinction are needed.
Geographic range size reduction is one of the main factors that increase extinction risk (Purvis et al., 2000;Saupe et al., 2015). Range contractions reduce the total population size, decreasing genetic variation. For example, a study of alpine chipmunks showed declining genetic diversity with reducing range size over the last century (Rubidge et al., 2012). Another study showed decreased fitness due to reduced genetic diversity and inbreeding in a declining population of southern dunlins (Blomqvist et al., 2010). These effects can be exacerbated by population fragmentation, which reduces gene flow even further (Gilpin, 1987). Geographic range contractions are suggested to fall into one of two models: the demographic hypothesis, which predicts that the range implodes toward the core of the historical range, and the contagion hypothesis, which expects the range to shrink away from local negative impacts (Channell & Lomolino, 2000) (Figure 1). The demographic hypothesis assumes that the intensity of environmental change is uniform across the range, whereas the contagion hypothesis assumes a single, local extinction factor, with chances of survival increasing further away from the source of the perturbation (Channell & Lomolino, 2000). It is not clear which of these hypotheses applies more commonly in nature. Range contractions preceding modern extinctions are often best described by the contagion hypothesis (e.g., Boakes et al., 2018;Lomolino & Channell, 1995;Rando et al., 2011;Steadman et al., 2005) but in the majority of these cases human influences were the cause of extinction. To understand the drivers of extinctions where humans are not the main agent of forcing, we need to study range size dynamics further back in time.
Declining populations often experience changes in morphology, either as a result of reduced genetic and/ or morphological variation, through inbreeding or otherwise, or as a direct response to external biotic or abiotic pressures. Body size has been particularly well studied in this regard, with reductions in average size reported at both the community and species levels. For example, average size in macroinvertebrate assemblages has been found to decrease immediately prior to the end-Permian mass extinction (Kiessling et al., 2018;Wu et al., 2018;Zhang et al., 2016Zhang et al., , 2017 and up to several hundred thousand years before the early Toarcian extinction event (García Joral et al., 2018;Piazza et al., 2019). Pre-extinction dwarfing of individual species has been noted in planktonic foraminifera (Wade & Olsson, 2009) and Pleistocene megafauna (Guthrie, 2003). Shifts in body size can precede population decline, which has recently led several studies to conclude that a combination of population abundance and body size is a more reliable early warning signal of population collapse than population size alone (Baruah et al., 2019;Clements & Ozgul, 2016). Baruah et al. (2019) found that this is most likely when environmental change was slow relative to generation time. Trait changes can buffer the negative effects of environmental change for a while, until a critical threshold is crossed and adaptation is no longer possible. Thus, when studied quantitatively, traits such as body size provide a promising means to study extinction dynamics in deeper time, especially because environmental change in the geological past tends to be slow relative to generation time (the immediate environmental consequences of meteorite impacts are the most obvious exception). Studies of geographic range, population abundance, and morphology through time and space require highly resolved fossil records. Most fossil records lack this resolution, but marine microfossils are exceptional. Their large population size, global distribution, and excellent preservation potential make it possible to generate abundance records and morphological reconstructions at resolutions exceeding environmental change over millions of years. Additionally, environmental reconstructions from the same sample medium allow for direct comparisons of biotic and abiotic change. Here, we study range size, abundance and morphology of the planktonic foraminifer Globoconella puncticulata preceding its extinction. This globally distributed Pliocene species became extinct around 2.4 Ma (Wei, 1994), shortly after Earth's last major shift in global climate state: the intensification of Northern Hemisphere glaciation (iNHG, ∼3.5-2.5 Ma) (Bailey et al., 2013;Westerhold et al., 2020). During this interval, Earth transitioned from a unipolar glacial state, with large ice sheets only on Antarctica (Spray et al., 2019), to one that was sufficiently cool by ∼2.5 Ma to induce growth of large ice sheets in the Northern Hemisphere (Shakun et al., 2016) and associated rafting of icebergs across the subpolar North Atlantic Ocean (Shackleton et al., 1984). Previous work on G. puncticulata has demonstrated diachronous population extinctions around the world (Scott et al., 2007), a decrease in test size prior to extinction in a mid-latitude population (Brombacher et al., 2017a) and cyclical variations in abundance corresponding to glacial-interglacial cycles in a subpolar population (Chapman et al., 1998). We studied 16,091 specimens of G. puncticulata from six Atlantic populations, over a 600,000-years time interval preceding its extinction at ∼2.4 Ma to investigate the consistency of these results along a latitudinal transect (Figure 2a). In particular, we address the following questions: 1. Does the geographic range contract before extinction?
We investigate extinction timing across populations. Where extinction is diachronous, we test the competing predictions of the demographic hypothesis (ranges implode toward the core of the historical range) and the contagion hypothesis (range shrinks away from local negative impacts).
2. Does shell size decrease before extinction?
One mid-latitude population of G. puncticulata showed a decrease in size prior to extinction (Brombacher et al., 2017a). Here, we investigate whether size changed consistently across populations to test whether size changes are a reliable early warning signal of extinction in other/all populations.
3. Do populations respond to the same environmental parameters similarly? BROMBACHER ET AL.   (Lisiecki & Raymo, 2005) with key Marine Isotope Stages (MIS) labeled, and ice rafted debris (IRD, orange) supplied to the open ocean by icebergs calved from ice sheets large enough to reach the coastlines, at Site U1313 (Bolton et al., 2010;Lang et al., 2014) in our study interval (b) and oxygen isotopes in the past 5 million years (Lisiecki & Raymo, 2005) with the temporal range of the study species (red) (c).
The demographic hypothesis assumes that response to environmental change is uniform across the species' range. Here, we compare local climate reconstructions to individual populations' abundance and size to test if this was the case, and to determine whether any differences among populations can be explained by differential response to local environmental perturbations.

Study Species
The planktonic foraminifer species G. puncticulata (Deshayes, 1832) first appeared in the early Pliocene around 4.6 Ma and became extinct at 2.41 Ma (Wei, 1994), shortly after the Neogene/Quaternary boundary (∼2.58 Ma) and the culmination of the late Pliocene-earliest Pleistocene intensification of Northern Hemisphere glaciation at ∼2.5 Ma (Bailey et al., 2013). The species was globally distributed with highest abundance in thermocline waters at middle and low latitude sites (Aze et al., 2011;Kennett & Srinivasan, 1983). A review by Scott et al. (2007) of first and last occurrences in Atlantic, Mediterranean, Pacific, and Southern Ocean sites shows asynchronous population extinctions of G. puncticulata, suggesting local environmental controls on population survival. Chapman et al. (1998) reported cyclical variations in relative abundance of a subpolar population in response to glacial/interglacial cycles, and a mid-latitude Atlantic population showed a decrease in test size ∼300,000 years prior to extinction (Brombacher et al., 2017a), possibly indicating pre-extinction dwarfing (Wade & Olsson, 2009).

Study Sites
We studied specimens of G. puncticulata from six Atlantic Ocean sites ( Figure 2) spanning ∼83° latitude: • Ocean Drilling Program (ODP) Site 981 (Jansen et al., 1996): A subpolar Northern Hemisphere site (55°N) currently situated under the North Atlantic Current. • Integrated Ocean Drilling Program (IODP) Site U1313 (Channell et al., 2006): A mid-latitude Northern Hemisphere site (41°N) situated at the northern edge of the oligotrophic North Atlantic Subtropical Gyre. • Deep Sea Drilling Program (DSDP) Site 606 (Ruddiman et al., 1987): A mid-latitude Northern Hemisphere site located just south of Site U1313 within the North Atlantic Subtropical Gyre (37°N). Site U1313 and Site 606 are situated closest to one another geographically and are the most similar of the sites studied oceanographically. • ODP Site 659 (Ruddiman et al., 1988): A subtropical Northern Hemisphere site (18°N) situated in the Canary Current off the Northwest African margin. • ODP Site 925 (Curry et al., 1995): A western equatorial site (3°N) located in the Southern Equatorial Current transporting South Atlantic surface waters fed from the Benguela Current into the North Atlantic. This site experiences high export production due to nutrient-rich Amazon river water runoff. • ODP Site 1264 (Zachos et al., 2004): A subtropical Southern Hemisphere site (28°S) located at the eastern edge of the South Atlantic Subtropical Gyre in the Benguela Current. Of all study sites, this site is least directly affected by climate feedbacks associated the intensification of Northern Hemisphere glaciation (Zachos et al., 2004).

Abundance and Geographic Range Size
All samples were dry-sieved over a >150 μm mesh sieve and divided using a Green Geological aluminum microsplitter until a single split contained 50-100 adult individuals of G. puncticulata, all of which were picked to avoid size biases. About 16,091 specimens of G. puncticulata were analyzed. To enable direct comparisons of abundance among sites, abundance per time slice was calculated as flux rather than relative abundance. The estimated total number of specimens in a sample (the number of individuals found in an analyzed split multiplied by the split fraction) was divided by sample weight, dry bulk density, and the total time contained in the sample with final numbers expressed as individuals/cm 2 /kyr.

Morphology
To measure test size, foraminiferal shells were mounted on glass slides in groups of 20 individuals using transparent double-sided adhesive tape with the apertures facing upwards. Groups were imaged using an Infinity 3 Lumenera camera mounted on an Olympus SZX10 light microscope, and test size was extracted from the images using an automated image analysis macro in the Image Pro Premier software. We measured size as test area. Previous error analyses of a subset of the data presented here have shown that test area is a repeatable measure (Brombacher et al., 2017b) and a more accurate representation of three-dimensional body size than test diameter (Brombacher et al., 2018a).

Response to Local Environmental Change
Abundance and morphology of one mid-latitude population of G. puncticulata have been shown to respond to a combination of environmental factors, rather than single variables (Brombacher et al., 2018b). Therefore, we study biotic response to multiple abiotic drivers. Planktonic foraminifera are known to have preferred temperature ranges (Lombard et al., 2009(Lombard et al., , 2011 and growth rates are influenced by food availability (Bé et al., 1981;Takagi et al., 2018). Additionally, ocean pH influences selection for larger shell size and thickness with decreasing pH (Hemleben et al., 1989).
To avoid offsets among local environmental records, it is important to compare population response to environmental change reconstructed using the same proxies. Two of our study sites met these criteria: Site U1313 has high-resolution reconstructions of sea surface temperature (SST) and export productivity available for our study interval, and SST and export productivity reconstructions from Site 982 were used as indicators of local environment at Site 981 given the close proximity of these two sites (Jansen et al., 1996).
Reconstructions of sea surface temperature are based on the alkenone unsaturation index U kʹ 37 (982: Lawrence et al., 2009;U1313: Naafs et al., 2010 and productivity records are reconstructed from C 37 alkenone mass accumulation rates (C 37 MAR) (982: Bolton et al., 2011;U1313: Naafs et al., 2010. Additionally, we used a record of atmospheric CO 2 (Martínez-Botí et al., 2015) to compare biotic response to global as well as local environmental change. Where climate data point ages were offset relative to our new biotic data, we followed the approach of Brombacher et al. (2018b) using Generalized Additive Models (GAMs) to interpolate the climate records to the foraminiferal sample ages by applying a smooth function to the climate data (R package "mgcv," version 1.8.33 (R Core Team, 2013;Wood, 2006)). The response variables in a GAM depend linearly on the additive combination of polynomial functions multiplied by statistically estimated coefficients (Wood, 2006) making them more suitable to interpolate environmental data (which often do not vary linearly with time) than simple linear interpolations. This approach is capable of distinguishing the glacial-interglacial dynamics of our study interval (Brombacher et al., 2018b). Although including only two of our six study sites in this analysis will not allow us to quantify the full range of responses of G. puncticulata to climate change across its geographical range, this is not necessary to distinguish between the two competing hypotheses of the demographic and contagion hypotheses. The two sub-Arctic and mid-latitude sites experienced the largest environmental change across the species' range due to their proximity to the Northern Hemisphere ice sheets (Herbert et al., 2010;Lawrence et al., 2009;Naafs et al., 2010), providing an opportunity to study biotic response to the maximum environmental change experienced by both a core and edge of range population. The multiproxy climate records at these two sites allow investigation of synergistic climate effects (Brombacher et al., 2018b) to provide a first indication of possible variation among populations.
Environmental conditions at Sites U1313 and 982 were influenced by both global and local factors. SST records at both sites show increasing glacial-interglacial variation through our study interval but, whereas Site 982 also experienced an overall decreasing trend in SST, only average glacial temperatures decreased at Site U1313 (Lawrence et al., 2009;Naafs et al., 2010). Changes in ocean circulation and global atmospheric CO 2 decline were responsible for SST fluctuations at both sites (Friedrich et al., 2013;Lawrence et al., 2009;Naafs et al., 2010), and glacial-interglacial SST differences at Site 982 were further amplified by the site's close proximity to the Northern Hemisphere ice sheets (Lawrence et al., 2009). Transient peaks in export productivity are caused by increased wind-driven mixing of deeper, nutrient-rich waters at Site 982 (Bolton et al., 2011) whereas southward deflections of the sub-Arctic front and increased glacial dust delivery enhanced productivity at Site U1313 (Bolton et al., 2018;Naafs et al., 2010Naafs et al., , 2012. Our study sites are situated in different paleoceanographic settings that respond differently to environmental change (Figure 2a). To test the similarities and differences in population response to the same underlying environmental driver, the intensification of Northern Hemisphere Glaciation (iNHG), and its local expression we used a Linear Mixed Effect model using the R package "nlme" (version 3.1.152) (Pinheiro et al., 2021;R Core Team, 2013). This model looks for linear trends among subsets of the data (here: populations) by comparing support for six different scenarios: 1. no ocean-wide climate trends, random effects, or variance differences between sites (model input: trait ∼ 1) 2. no ocean-wide climate trends or random effects but variance differences between sites (model input: trait ∼ 1|site) 3. no ocean-wide climate trends, but site-specific intercept differences and variance differences (model input: trait ∼ 1 + (1|site)) 4. no ocean-wide climate trends, but site-specific intercept differences, temperature relationships, and variance differences (trait ∼ 1 + (sst|site)) 5. ocean-wide climate trends, site-specific intercept differences and variance differences (trait ∼ climate + (1|site)) 6. ocean-wide climate trends, site-specific intercept differences, temperature relationships, and variance differences (trait ∼ climate + (sst|site)) (see also Figure S1) Models without random effects were fitted using generalized least squares models. Variance differences among sites were assessed using the varIdent function in nlme (Pinheiro & Bates, 2000). The model with the lowest Akaike Information Criterion (AIC) value (a measure for model parsimony) and the highest Akaike weight (a measure for relative model likelihood) has the most statistical support among those models fitted. "Climate" was taken as the combination of local SST and export productivity and global atmospheric CO 2 , and their two-way interactions (model input: climate = (SST + productivity + CO 2 ) 2 , with the squared term calculating two-way interactions between the environmental parameters). A visual inspection of the residuals shows that they are approximately normally distributed, satisfying the model prerequisites ( Figures S2-S5).

Results
Abundance varies among sites (Figures 3a-3d The two populations at sites with environmental reconstructions available (Sites U1313 and 981) show different responses to local abiotic change ( Figure 5). The final extinction as well as two temporary extirpations at Site 981 occur during SST minima associated with MIS 96 (extinction), 100, and G6 (extirpations), whereas abundance at Site U1313 shows no response to local temperature or productivity. In contrast, mean test size at Site U1313 decreases during the SST minimum and productivity maximum associated with MIS G6 but no direct size response to temperature or productivity is observed at Site 981. The Linear Mixed Effects model results indicate no trend in the pooled abundance data of Sites U1313 and 981, but response to local environment differs between the subpolar and mid-latitude populations (    Size shows both a detectable trend in the pooled data, as well as site-specific trends of different magnitude (Scenario 6, Table 3). Unlike the abundance models, there is a smaller but still statistically significant impact of variance differences between the sites (cf. Models 1 and 2). The impact of the variance differences is, however, much smaller than including global trends in the climate variables and their interactions (cf. Models 3, 5, and 6)-that the biggest jump in AIC scores is through the climate variables as the fixed rather than random effects implies much stronger global controls on size than abundance (cf., Tables 2 and 3). There are clear site-specific differences in how these planktonic foraminifera respond to sea surface temperature changes (cf. the improvement from Model 5 to 6 in Table 3), supporting conclusions that the organismal response to local environmental settings is not constant through space and time.

Range Contraction
The range contraction of G. puncticulata preceding extinction does not clearly follow predictions consistent with either the demographic or the contagion hypothesis. The demographic hypothesis predicts a range contraction to the core (e.g., edge Sites 981 and 925 first, followed by core populations U1313, 606, 659, 1264), and while one core population (Site 659) is indeed the last to become extinct, the core populations at Sites U1313 and 606 disappear much earlier (Figures 4 and S6). Similarly, the contagion hypothesis predicts extinctions away from the presumed extinction driver; in this case, the environmental response to the expansion of Northern Hemisphere ice sheets (e.g., sub-Arctic Site 981 first, followed by mid-latitudes Sites U1313 and 606, subtropical Sites 659 and 1264, and finally tropical Site 925), and while the first two extinctions do follow this pattern it breaks down with the subsequent extinction of the equatorial population at Site 925. A combination of the two hypotheses better explains the data: the subpolar (Site 981) and equatorial (Site 925) edge populations disappear early on as predicted by the demographic hypothesis, but the extinction of BROMBACHER ET AL.  Note. ΔAIC represents the difference between AIC and the set's minimum AIC. The best performing model based on lowest AIC and highest Akaike weight is indicated in bold. See Figures S4 and S5 for residual distributions.

Table 3 Linear Mixed Effect Model Results of Globoconella puncticulata Size Response to Local Environmental Change
Northern Hemisphere core populations U1313, 606, and 659 follows a North-South pattern resembling the contagion hypothesis. These results suggest that population dynamics in core populations differ from those at the edge, and that biogeographic variation can obscure the detection of extinction driving mechanisms.
Our results differ from many studies of population range contractions. A review by Channell and Lomolino (2000) found that most modern ranges contract according to the contagion hypothesis. However, in most of these cases, the extinction pressure came from humans. For example, Steadman et al. (2005) showed that, for late Quaternary sloth populations, extinction matched human arrival rather than global climate change and extinction of modern Galliform populations occurred close to the core in a human-dominated landscape (Boakes et al., 2018). Boakes et al. (2018) also showed that, in the absence of humans, extinction was equally likely to occur near the core or the edge of the geographical range, offering support for both the demographic and extinction hypotheses. Our data strongly suggest that contagion and demographics both played a role in extinction, but data from more extensive geological intervals and different types of species are needed to test the generality of this result.
A combination of the demographic and contagion hypothesis, rather than either of the theoretical "endmember" extinction hypotheses on their own, might offer a more plausible explanation of extinction in general. The demographic hypothesis is based on the observation that populations are larger and more variable in the center of their range than at the edge, which makes edge populations the most vulnerable to extinction (Channell & Lomolino, 2000). To test if an extinction follows the demographic hypothesis, the extinction factor must be synchronously and uniformly distributed across space (Channell & Lomolino, 2000). However, except in rare cases of species with an extremely small range, environmental change is never homogenous across the range. If certain areas are affected more than others, these populations would be more likely to follow the contagion hypothesis, calling into question the extent to which extinctions will ever solely resemble the demographic hypothesis. Similarly, the contagion hypothesis predicts extinction closest to the external stressor, but core populations might be more resilient to external perturbations than edge populations, decreasing detectability while also introducing elements of the demographic hypothesis.

Abundance
In contrast to Chapman et al. (1998), we did not find clear cyclical variations in abundance as a response to glacial-interglacial cycles. This could be explained by different measures of abundance: we report abundance as absolute abundance or flux (individuals per area per kyr) whereas Chapman et al. (1998) measured abundance of G. puncticulata relative to other species in the assemblage. Relative abundance is a measure susceptible to closed sum problems. In the subpolar North Atlantic, decreases in the relative abundance of G. puncticulata during the early Northern Hemisphere glacials may well have been driven by increases in the abundance of cold-adapted species such as Neogloboquadrina pachyderma, Neogloboquadrina atlantica, and Neogloboquadrina incompta. Flux on the other hand is more susceptible to age model uncertainty than relative abundance. Based on a maximum uncertainty of ∼6 kyr in the LR04 age model in our time interval (Lisiecki & Raymo, 2005), and assuming one tie point per glacial-interglacial transition (i.e., one tie point per 20 kyr, see (Bolton et al., 2010)), the time contained in a sample and thus the flux could vary by a factor of ∼0.25. However, as our flux patterns vary by 2-3 orders of magnitude (Figures 3a-3d) the uncertainty introduced by the age model is unlikely to explain the absence of glacial-interglacial cyclicity in our data.
The use of relative abundance could also explain the different extinction timing we found for Site 659 as compared to Chapman et al. (1998). Chapman et al. (1998) reported extinction of the subtropical population at 2.44 Ma, similar to other populations, whereas we found G. puncticulata there for another 100,000 years until 2.34 Ma despite similar sampling resolutions. Assemblage counts typically analyze 300-500 specimens per sample, whereas here we included 50-100 specimens of the study species alone, or all specimens in a sample if there were <50. This approach is more likely to find (temporarily) rare species. Abundance of G. puncticulata decreased by a factor 10 on average in the last 100,000 years before extinction at Site 659, so it is possible that an assemblage study would have been less likely to find the species on the approach to extinction.

Pre-extinction Dwarfing
Two of the six studied populations (Sites U1313 and 606) show a sudden size decrease consistent with pre-extinction dwarfing as described by Wade and Olsson (2009) (Figure 3f). Specimens showing pre-extinction dwarfing in that study mostly came from core populations, so it could be that dwarfing does not occur uniformly across space. Another possibility is that individuals at edge sites are smaller than those at core sites, resulting in size changes that are more difficult to detect due to the relative importance of different climate drivers (Brombacher et al., 2018b). Our Linear Mixed Effect model results support this explanation, showing different responses of size to local environmental change between populations at Sites U1313 and 981. More data on synchronicity of size changes across populations are needed before firm conclusions can be drawn, but our results suggest that the expression of early warning signals may vary across populations. In particular, a lack of change in one population does not necessarily imply low extinction risk or an absence of early warning signals for the species as a whole given the importance of the local paleoceanographic context.

Causes of Extinction
The differences in population dynamics through space and time that we document suggest that the extinction of G. puncticulata cannot readily be explained by a single driver. If temperature was the main driver we would expect to see cyclical abundance variations in all populations with local extirpations during severe glacial stages. The subpolar population at Site 981 followed this pattern with temporary extirpations during glacials MIS G6 (2.72 Ma) and 100 (2.54 Ma, see Figure 5), but no other populations show similar responses. As predicted by both the demographic and contagion hypotheses, the population at Site 981 disappears first, but after that there is also no clear chain of extinction from North to South. Core Sites U1313, 606, and 659 disappear in latitudinal sequence, but the equatorial edge population (Site 925), which under the contagion hypothesis should thrive under cooler conditions resembling the core environment, became extinct shortly after the first mid-latitude population disappeared. These results suggest that individual populations develop local adaptations that do not translate simply across space.
Local adaptation is further supported by the results of our Linear Mixed Effect models, which show significant differences in population response to local environment. As we analyzed response to climate change in only two populations the species' full range of response is not quantified, but the different responses observed at two relatively close sites suggests that differential response among populations to the same environmental parameters is not uncommon. These results agree with previous studies that found individualistic population response to climate change (Levin, 2019;Stewart, 2009;Stewart et al., 2010). Stewart (2009) found that Quaternary populations respond individualistically to biotic interactions through time and space. Stewart et al. (2010) also showed that the location of refugia following range contraction depends on population-specific adaptations and is therefore not always located in the core of the initial range. Levin (2019) argued that existing range contraction hypotheses oversimplify extinction precisely because they do not take local adaptation into account. Instead, Levin (2019) proposes that distinct intraspecific lineages are the focal point of species' extinction, and that only by combining all intraspecific lineage responses can we understand the extinction of species.
Population-specific response to environmental change increases chances of survival of at least some populations (Levin, 2019). In our data, however, low abundance and suboptimal adaptation to local environmental change in edge populations likely led to an increasingly fragmented species distribution. Temporary extirpations of the equatorial population (Site 925) at 2.8 and 2.7 Ma (Figure 3d) resulted in an isolated Northern Hemisphere population with reduced global gene flow. The gradual extinction chain of the core populations briefly resulted in a subtropical refugia at Site 659 which persisted for another 30,000 years after the other populations had become extinct. The resistance to extinction at this site may be attributed to a combination of high abundance and phenotypic variability typical for a core population which reduces extinction risk (i.e., demographic hypothesis). Alternatively, being located furthest away from Northern Hemisphere ice influences of all Northern Hemisphere core populations (i.e., contagion hypothesis), or an optimal, site-specific microclimate may have made Site 659 a refuge for G. puncticulata. Populations in refugia, no matter how well adapted initially, are vulnerable to extinction because of their small population size and low genetic diversity, and absence of gene flow and migration from other areas following population decline (Blomqvist et al., 2010;Rubidge et al., 2012). Abundance at Site 659 is highly variable in its final 100,000 years (Figure 3c), suggesting increased population instability before extinction.

Directions for Future Work
To understand extinction dynamics in deep time, data from many more microfossil groups are needed. Studies of extinction through time and space require high-resolution fossil records from multiple locations with high abundance of complete fossils. Most fossil records lack this resolution, but those of marine microfossils are an exception. The application of marine microfossils to biostratigraphy has already provided insights in spatial dynamics preceding extinction: it is well known that some species go extinct almost simultaneously around the world and can be used as global biohorizons (e.g., Falzoni et al., 2018;Gibbs et al., 2005;Huber et al., 2017;Wade et al., 2011) whereas others undergo extinction diachronously and are of local stratigraphic use only (e.g., Du et al., 2020;Falzoni et al., 2018;Gibbs et al., 2005;Huber et al., 2017;Kawagata et al., 2005;Prell & Damuth, 1978). However, these patterns are rarely extensively studied from a biological perspective. Global climate change is occasionally mentioned as a cause of varying population extinction timings (Gibbs et al., 2005;Prell & Damuth, 1978), but range contractions and individual population response to local environmental change are rarely quantified. Here, we have provided an example of how data on microfossils and climate can provide insights into population dynamics preceding extinction in geologic time.

Conclusion
The extinction of G. puncticulata was preceded by changes in the geographic range, abundance and body size that varied among populations. Range contraction prior to extinction showed characteristics of both the demographic (contraction toward the core) and contagion hypotheses (shrinking away from unfavorable conditions). The two northern-most populations disappeared first in keeping with the contagion hypothesis but a subtropical core population outlasted a tropical edge population before ultimately becoming extinct, suggesting influences of the demographic hypothesis as well. These two broad competing hypotheses are very likely a false dichotomy in many cases-here, the response to local environmental change varied between a subpolar and a mid-latitude population, supporting previous inferences that populations are differentially adapted to temperature, nutrient availability and other microclimatic factors. Ultimately, global extinction can only be understood by unraveling the dynamics of individual populations and by understanding how their dynamics are linked through space and time. To test the generality of our results, more data are needed from foraminifera and other fossil groups through different time intervals and climatic settings. The method of co-registered marine microfossil and paleoclimate data is well-poised to yield transformative new insights into extinction dynamics.