Evaluating the responses of three closely related small mammal lineages to climate change across the Paleocene–Eocene thermal maximum

Abstract

Interpreting the impact of climate change on vertebrates in the fossil record can be complicated by the effects of potential biotic drivers on morphological patterns observed in taxa. One promising area where this impact can be assessed is a high-resolution terrestrial record from the Bighorn Basin, Wyoming, that corresponds to the Paleocene–Eocene thermal maximum (PETM), a geologically rapid (~170 kyr) interval of sustained temperature and aridity shifts about 56 Ma. The PETM has been extensively studied, but different lines of research have not yet been brought together to compare the timing of shifts in abiotic drivers that include temperature and aridity proxies and those of biotic drivers, measured through changes in floral and faunal assemblages, to the timing of morphological change within mammalian species lineages. We used a suite of morphometric tools to document morphological changes in molar crown morphology of three lineages of stem erinaceid eulipotyphlans. We then compared the timing of morphological change to that of both abiotic and other biotic records through the PETM. In all three species lineages, we failed to recover any significant changes in tooth crown shape or size within the PETM. These results contrast with those documented previously for lineages of medium-sized mammals, which show significant dwarfing within the PETM. Our results suggest that biotic drivers such as shifts in community composition may have also played an important role in shaping species-level patterns during this dynamic interval in Earth history.

Introduction

A primary challenge to understanding the evolutionary consequences of climate change is separating the effects of diverse proximate drivers, such as changing biotic interactions between competitors and abiotic water availability, on individual species (Cahill et al. 2013). Sampling biases in the fossil record combined with lack of direct control of variables, as one would have in a standard laboratory experiment, makes separating the signal of abiotic from biotic drivers difficult, if not impossible (Gill et al. 2009; van Gils et al. 2016). In certain records, one possible solution is to evaluate the synchronicity of patterns associated with distinct hypothesized drivers. If changes in the fossil record and patterns of different drivers are found to be asynchronous, then the relative timing of change in a species can potentially be used to differentiate abiotic from biotic drivers (Gill et al. 2009). Examples of such fossil records are rare (Gill et al. 2009; Terry 2018). One such example comes from a sustained collaborative research effort to document high-resolution floral and faunal biostratigraphy correlated with similarly sampled chemostratigraphy across and through the Paleocene–Eocene thermal maximum (PETM) within the southeastern Bighorn Basin (BHB) of Wyoming (e.g., Strait 2001; Wing et al. 2005; Yans et al. 2006; Kraus and Riggins 2007; Smith et al. 2007; Chester et al. 2010; Adams et al. 2011; Rose et al. 2011, 2012; Secord et al. 2012; Baczynski et al. 2013, 2017; Kraus et al. 2013; Bourque et al. 2015; Morse et al. 2019). Abiotic and biotic records show the kind of asynchronous pattern that allows for their differentiation during the PETM. Within the boundaries of this event, abiotic variables such as temperature and hydrology continued to change locally (Secord et al. 2012; Kraus et al. 2013; Baczynski et al. 2019). In contrast, the fossil record of flora and fauna during the PETM preserves relatively static species composition, with turnover largely confined to the PETM boundaries (Koch et al. 1992; Gingerich 2001; Wing et al. 2005). We use this fossil record to document morphological change in ecologically similar lineages of small-bodied mammals, classified as three distinct lineages of closely related stem erinaceid (“erinaceomorph”) eulipotyphlans, as a test case to address whether abiotic climate change is the primary driver that explains morphological changes in mammalian lineages across and through the major global climate change of the PETM.

Categorizing Drivers and Responses

Studies documenting how species are influenced by climate change have often, reasonably, focused on measuring direct responses to abiotic climate (e.g., Kurtén 1960; Smith et al. 1995; McGuire 2010; Prothero et al. 2012; Williams and Blois 2018). However, climate change affects many aspects of an ecosystem simultaneously (Davis et al. 1998; Alexander et al. 2015; Eskelinen et al. 2017; Mokany et al. 2019). A change in climate may instigate additional, interdependent drivers, both abiotic and biotic, that can directly influence the survival of a species (Davis et al. 1998; Cahill et al. 2013). This range of potential proximate drivers introduces complexity into predictions of how any single species will ultimately respond to climate change (Graham 1986; Roy and Pandolfi 2005; Stewart 2008). More recent studies have provided new research frameworks to help tease apart the variety of factors driving responses to climate change, both in the modern biota and the fossil record (Gill et al. 2009; Cahill et al. 2013; van Gils et al. 2016).

Proximate drivers can be broadly categorized as abiotic and biotic (Cahill et al. 2013). Abiotic drivers refer to climatic conditions, such as temperature, precipitation, and humidity, as well as factors such as salinity or acidity in aquatic or soil environments (Cahill et al. 2013). These factors can drive responses through physiological limitations of species. For example, lethal temperatures are related to body mass in woodrats (Neotoma), especially in arid climates (Smith and Betancourt 2006), and can explain shifts in body mass corresponding to temperature change during the late Pleistocene and Holocene of the Great Basin and Colorado Plateau (Smith et al. 1995).

Biotic drivers refer to interactions with other species. For heterotrophs, these interactions can include food types and abundance, competitors, predators, and species that contribute to forming microhabitats (Cahill et al. 2013). These factors can drive responses by forcing ecological limits on species survival. For example, the isotopic evidence of stronger Holocene niche partitioning in small mammals of the Smoke Creek Desert is better explained by changes in the resource pool, particularly the spread of invasive cheatgrass, than changes in paleoclimate (Terry 2018). Biotic factors are a frequent cause of local population declines that are ultimately related to climate change (Cahill et al. 2013; Westover and Smith 2020).

Categorizing drivers as abiotic versus biotic does not render them mutually exclusive when it comes to their role in influencing evolution or stasis in any single species. As experimentally documented, responses to warming can be the result of interactions between both physiological limitations and species interactions (Alexander et al. 2015; Eskelinen et al. 2017). Both categories of drivers can ultimately be driven by change in abiotic climate (Cahill et al. 2013). The utility of distinguishing the two lies in the ability to evaluate whether one or both categories are necessary to explain an apparent species’ response.

The PETM within the BHB

The terrestrial record of the PETM in the BHB, Wyoming, is one case in which categories of proximate drivers can be discriminated (Gingerich 2006; McInerney and Wing 2011). The PETM occurred at the beginning of the Eocene epoch, ~56 Ma (Gradstein et al. 2012). The onset of the PETM is recognized by a rapid negative carbon isotope excursion, which marks the beginning of the Eocene and indicates an increase in atmospheric carbon dioxide associated with large-scale range shifts of floras, as well as a fundamental reorganization of the mammalian biota of North America (Clyde and Gingerich 1998; Gingerich 2006; Woodburne et al. 2009; McInerney and Wing 2011; Wing and Currano 2013).

Globally, the climate of the PETM consists of a ~13,000 year warming “onset” interval followed by a ~115,000 year “body” period of relatively stable climate, then a ~42,000 year cooling “recovery” interval (Bowen et al. 2006; Aziz et al. 2008; Westerhold et al. 2009; McInerney and Wing 2011). Over the course of the PETM, global mean annual temperature changed by ~5–8°C. Locally to the BHB, warming continued from the onset into the early part of the body of the PETM, and cooling is well constrained within the recovery (Secord et al. 2012; Fig. 1A).

Figure 1. Patterns of M₁ crown area through the Paleocene–Eocene thermal maximum (PETM) (A) in comparison to the oxygen isotope record reconstructing mean annual temperature from δ¹⁸O of tooth enamel in the mammal Coryphodon. B, Patterns of M₁ crown area in the equid lineage Sifrhippus-Arenahippus (B) and in three ecologically similar erinaceomorphs (C–E): Colpocherus (C), Macrocranion (D), and Talpavoides (E). Warmer temperatures in A correspond to more positive δ¹⁸O values. A, B, Data from Secord et al. (2012). Dashed lines indicate the beginning and end of the carbon isotope excursion delimiting the PETM. Point shapes and colors indicate sampling bins. Brackets indicate comparisons. Asterisks next to brackets indicate statistically significant differences. (Color online.)

The PETM is also marked by major paleohydrologic transitions in the intermontane basins formed during the Laramide orogeny, including the BHB (Foreman et al. 2012; Kraus et al. 2013; Baczynski et al. 2017). The timing of paleohydrologic changes largely follows the timing of paleotemperature changes (Secord et al. 2012; Kraus et al. 2013; Baczynski et al. 2019). During the PETM onset and into the early PETM body, climate in the BHB became warmer and more arid (Kraus et al. 2013). It remained arid with more seasonal or episodic rainfall through the mid-body of the PETM (Wing et al. 2005; Foreman et al. 2012; Kraus et al. 2013; Baczynski et al. 2017, 2019). In the PETM recovery, a wetter climatic regime began, eventually returning the basin to climatic conditions approximating those of the late Paleocene.

The terrestrial record of the PETM in the BHB overcomes many of the practical limitations that usually prevent discrimination of abiotic and biotic drivers of change in the fossil record (Barnosky 2001). Previous work documenting the isotope- and paleosol-based climatic reconstructions provide a local abiotic record that can be compared with the biotic fossil plant and vertebrate records (Wing et al. 2005; Secord et al. 2012; Baczynski et al. 2013, 2017; Kraus et al. 2013). A high-resolution stratigraphic framework (Baczynski et al. 2013) coupled with large-scale screen-washing efforts has resulted in the recovery of many teeth of small-bodied vertebrate species from multiple time intervals before, during, and after the PETM across 16 km in the southern BHB (Baczynski et al. 2013). In several cases, these efforts produced sample sizes that are large enough to quantify simultaneous responses in multiple taxa (Polly 2003; Secord et al. 2012). Each species is also associated with functional traits that can help identify functional groups of species with broad similarities to one another in terms of their probable physiological tolerances and ecological interactions, and therefore are predicted to respond similarly to climate change drivers in general (Reed et al. 2006; Eronen et al. 2010; Terry et al. 2011). Finally, all records, both abiotic and biotic, are tied into the same stratigraphic framework as the vertebrate fossil record, allowing for a comparison of relative timing between climatic, floral, and faunal records (Baczynski et al. 2013).

One major advantage of the PETM record in the BHB is the degree to which the fauna has been studied, including mammals, terrestrial invertebrates, and reptiles (e.g., Holroyd et al. 2001; Gingerich 2006; Smith 2009; Smith et al. 2009; Wing and Currano 2013; Bourque et al. 2015). The mammalian fauna in particular has received multiple monographic treatments in addition to taxon-specific studies, forming a strong foundation for alpha taxonomy and species identifications (e.g., Gingerich 1989; Strait 2001; Smith et al. 2002; Heinrich et al. 2008; Chester et al. 2010; Rose et al. 2012). Biostratigraphic and biogeographic studies have also been undertaken (e.g., Gingerich 2001; Bowen et al. 2002; Gingerich and Smith 2006; Smith et al. 2006; Chew 2009; Solé and Smith 2013; Rankin et al. 2015; Morse et al. 2019; Fraser and Lyons 2020; van der Meulen et al. 2020), allowing for assessments of the nature and timing of mammalian community turnover. Floral and faunal turnover generally occurs at the biostratigraphic boundaries bracketing the PETM, but turnover is low or absent within the PETM itself as abiotic climate continued to change (Koch et al. 1992; Gingerich 2001; Wing et al. 2005; although see Hooker 2015). At the onset of the PETM, rapid floral turnover is associated with increased diversity and amount of insect damage, indicating a potential change in abundance of food for insectivorous mammals (Currano et al. 2016). The plant and insect communities turned over again through the PETM recovery (Wing and Currano 2013). By the end of the PETM, incidence of insect damage was similar to that of pre-PETM levels (Currano et al. 2016). Mammalian faunal change largely occurs either within the first 4–5 m above the PETM onset or at the biostratigraphic boundary that occurs mid-recovery (Bowen et al. 2001; Gingerich 2001; Gingerich and Smith 2006). Most importantly for this study, within the onset and body of the PETM of the BHB, documented community composition was relatively unchanging. In addition, post-PETM faunal composition was similarly stable for the next ~500 kyr, with relatively static levels of diversity and relatively few first or last appearances (Chew 2009). The biotic record in the PETM is reminiscent of a nonlinear community response to climate change, wherein community composition shifted at the beginning and end of the PETM but did not continue to change dramatically or in sync with further climatic change within the PETM (Lavergne et al. 2010; Walther 2010; Willis et al. 2010), although further direct study is required to assess this idea.

Although floral and faunal turnover is concentrated at PETM boundaries, it is unclear whether other kinds of responses, particularly changes in morphology potentially tied to functional shifts, follow similar patterns. Previous studies of similar questions have not had the temporal resolution to address the possibility of changes within the PETM itself (Clyde and Gingerich 1994; Wood et al. 2007). Within the context of the PETM, the temporal asynchrony between abiotic and biotic change documented thus far results in three expectations: (1) Biotically driven responses to climate change should be limited to the boundaries of the PETM. Specifically, a species-level morphological response should occur during the transition from the pre-PETM to early PETM, and from the late PETM to post-PETM. (2) Biotically driven responses should not occur within the PETM, when biotic community composition was relatively unchanging. (3) Abiotically driven responses may occur both within the PETM and across PETM boundaries, as abiotic drivers such as temperature and precipitation continued to change (Fig. 1A).

For the purposes of hypothesis testing, this framework makes observations of change in species characteristics across the boundaries of the PETM interesting, but not useful for distinguishing abiotic and biotic drivers of morphological change. Observations of significant change in species morphology across boundaries could be consistent with either category of driver, because both abiotic and biotic conditions changed across each boundary. However, observing species morphology across the boundaries is still important. For example, morphological stasis across the PETM boundary is inconsistent with a response to abiotic drivers, because temperature and hydrology changed markedly across the PETM boundary (Foreman et al. 2012; Secord et al. 2012; Kraus et al. 2013; Fig. 1A). Stasis across the PETM boundary could be due to a biotic driver if the most important biotic driver for a given species also did not change between our latest pre-PETM and earliest PETM samples.

This Study

Here, we study three species of small-bodied mammals as a test case for using asynchrony to address the overall hypothesis that abiotic climate change can explain mammalian responses to the PETM. The biological research framework of the fundamental niche, or the abiotic environmental factors permitting species persistence, and its applications, including ecological niche modeling (Soberón and Peterson 2005) and the use of vertebrate morphology as paleoclimate proxies (Dayan et al. 1991; Head et al. 2009; McGuire 2010), as well as results from previous research on the topic (Smith et al. 1995; Secord et al. 2012), predict that such a direct relationship between abiotic climate change and species morphological responses exists (Polly et al. 2011; Lawing et al. 2017).

We focus on small-bodied mammals because they are thought to be especially sensitive to climate change due to their narrower physiological tolerances (Barnosky et al. 2004; Terry and Rowe 2015). The narrower tolerances and sensitivity may make them likely to have a direct response to abiotic climate change, though examples in certain extant systems could also support expectations of a biotic response (Yom-Tov and Yom-Tov 2005; White and Searle 2007). We chose species from among the most common small-bodied mammals recovered within our study area in order to amass as much statistical power as reasonably possible and avoid conflating lack of power to detect change with an inference of lack of change. The three species studied here are stem members of the Erinaceidae (also known as Erinaceomorpha), extinct relatives of hedgehogs and moon rats (Novacek et al. 1985; Gunnell et al. 2008; Beard and Dawson 2009; He et al. 2012). Postcranial fossils from the BHB support the interpretation that the three studied species were terrestrial or scansorial (Penkrot et al. 2008; Manz et al. 2015; Penkrot and Zack 2016). Eulipotyphlan morphology can track changes in temperature, precipitation, and food availability on a decadal scale (Yom-Tov and Yom-Tov 2005; Poroshin et al. 2010). Therefore, there should be no detectable lag between environmental change and potential morphological responses in the fossil record based on this capacity for rapid change.

We measure the response to drivers in terms of relative shifts in size and shape of the lower molars of each species. Molars can be used to measure valuable proxy data for inferring body size and diet (Kay 1975; Gingerich et al. 1982; Legendre 1986; Hlusko et al. 2006; Pineda-Munoz et al. 2017; Thiery et al. 2017) that are fundamental aspects of species ecology and are likely to be perturbed by climate change (Bergmann 1848; Schmidt-Nielsen 1984; Cahill et al. 2013). Body size fundamentally constrains physiology and both biotic and abiotic interactions (Calder 1983; Schmidt-Nielsen 1984). Rapid shifts in body size that correspond with intervals of climate change have been empirically documented (Smith et al. 1995; Secord et al. 2012; Teplitsky and Millien 2014; van Gils et al. 2016). Dietary resources are the primary way that organisms interact with the environment in terms of their dental morphology (Kay 1975; Sheine and Kay 1977; Gingerich et al. 1982; Coiner-Collier et al. 2016). Food availability is an important biotic cause of species responses to climate change (Cahill et al. 2013; van Gils et al. 2016). Conversely, changes in diet may result from a physiological response to abiotic climate that affects body size and foraging behavior (Pyke 1984; Dickman 1988). Therefore, although we cannot measure the phenome of each species to document all possible morphological responses, a focus on molar size and shape captures relevant information about some of the most important systems linking an animal to its changing environment.

One challenge in characterizing response through morphological change is evaluating whether absence of evidence can be considered evidence of absence (Levinton 1983; Barnosky 1993). Even with an explicitly limited focus on a single morphological system such as molar crown surface morphology, amassing both the statistical and methodological power to avoid overlooking potential change is a nontrivial challenge (Wood et al. 2007). To overcome this challenge and exhaustively explore potential responses in terms of molar crown surface morphology, we illustrate a workflow for the study of shape that leverages the utility of three-dimensional, univariate, angular, and topography morphometric tools to study complex morphology. In characterizing complex morphology, there is a trade-off between sample size and the dimensionality of our measurements. Three-dimensional measures of morphology that characterize the entire tooth are desirable, because they capture whole-crown morphology without being limited to the specific aspects of crown shape captured by univariate measurements (Boyer et al. 2015). They allow for the possibility of discovering unanticipated patterns of shape change (Zelditch et al. 2012). However, often only a limited number of specimens meet the stringent quality thresholds of being unworn and having the completely intact crowns necessary for meaningful 3D measurements and consistent comparisons among individuals (Vitek et al. 2017). In contrast, univariate linear and angular measurements can be taken from a much larger sample of partial and worn specimens, but the choice of where to focus linear measurement efforts may result in failure to capture important aspects of shape (Adams et al. 2004). Our multistep approach to quantifying shape leverages the detail of high-throughput 3D morphometric tools as well as the larger sample sizes of lower-dimensional tools. We use 3D analyses to inform univariate measurement choices, then use results from both types of analyses as complementary characterizations of shape through time. As part of characterizing morphology, we integrate function through analyses of shape along the wear surfaces that form in the process of chewing as well as metrics of whole-tooth topography that correlate with dietary variation (Kay and Hiiemae 1974; Ungar et al. 2018).

Materials and Methods

Studied Specimens

To identify the most common species of small mammals in the PETM sections, we comprehensively reviewed well-sampled screenwash collections from before, during, and after the climate event from the BHB. Out of a dataset of more than 21,000 vertebrate fossils collected from a single stratigraphic section in the BHB that spans the latest Paleocene to the early Eocene interval, including the PETM, we were able to identify 980 lower cheek teeth (P₄–M₃) of eulipotyphlans, the most commonly represented clade (mostly isolated teeth). Both left and right isolated teeth were sampled, because alluvial processes fragment and remove so much of the skeleton when depositing small mammal fossils such as these (Korth 1979) that we felt the added gain in sample size outweighed the relatively low risk of sampling the right and left molars of the same individual. Institutional abbreviations for specimens are as follows: UF, Florida Museum of Natural History, Vertebrate Paleontology, Gainesville, Fla.; USNM, National Museum of Natural History, Washington, D.C.

To classify morphotypes, we used the literature, comparative specimens, well-preserved specimens from the interval containing multiple tooth positions, and an exploratory morphometric approach to develop suites of diagnostic characters (for example workflow, see Vitek et al. 2017). Among the 980 cheek teeth, all were classified into distinct groups without any intermediates, indicating that these methods were robust at diagnosing morphotypes. The three most common types of stem erinaceids in Wa-0 were post hoc identified as the following species: cf. Colpocherus sp. (hereafter, Colpocherus), Macrocranion junnei (hereafter, Macrocranion), and Talpavoides dartoni (hereafter, Talpavoides; Bown and Schankler 1982; Smith et al. 2002; Rose et al. 2012). Further detail about the morphotyping process and its results is being written as part of a companion study describing the taxonomy and systematics of stem erinaceids within the section.

The studied time interval contains multiple mammalian biozones, representing one primary way of binning time into ecologically coherent units. All biozones are recognized by first and last appearances of larger-bodied species (none of them eulipotyphlans), eliminating circularity between relative divisions of time and identification of the focal species in this study. The latest Paleocene pre-PETM interval is contained within the Clarkforkian Cf-3 biozone. It is marked by the first occurrence of the phenacodontid condylarth Copecion and ends with the first appearance of the phenacodontid condylarth Meniscotherium (Gingerich 2001; Secord et al. 2006). The beginning of the PETM is coincident with the start of the Wasatchian North American Land Mammal Age (Koch et al. 1992) and the Eocene Epoch (Aubry et al. 2007). The earliest Wasatchian biozone, Wa-M, is marked by the first occurrence of Meniscotherium priscum and contains the first 4–5 m of the PETM onset (Bowen et al. 2001; Gingerich 2001; Gingerich and Smith 2006). It is difficult to recognize in the southern BHB, where the first occurrence of Wasatchian mammals frequently includes those characteristic of the next oldest biozone, Wa-0 (Rose et al. 2012; Baczynski et al. 2013). The Wa-0 biozone is recognized by the first appearance of the equid Sifrhippus sandrae and continues through much of the rest of the PETM, including the remainder of the onset, the entire body, and the beginning of the recovery intervals (Bowen et al. 2001; Gingerich 1989, 2001). The end of the PETM recovery shortly follows the first appearance of the perissodactyl Cardiolophus radinskyi, which defines the beginning of the following biozone, Wa-1 (Gingerich 2001). The Wa-1 biozone is preserved in only about 14 m of section in the southern BHB (Morse et al. 2019). To more fully characterize the post-PETM interval, we also include specimens from the subsequent Wa-2 biozone (Gingerich 2001; Table 1).

Table 1. Bins used in analyses in relation to geologic time, biostratigraphy, and changes in abiotic climate. MAT, mean annual temperature; NALMA, North American Land Mammal Age; PETM, Paleocene−Eocene thermal maximum. 1, Secord et al. 2006; 2, Secord et al. 2012; 3, van der Meulen et al. 2020; 4, Chew 2009. Note that van der Meulen et al. 2020 measured a shorter total duration of the PETM than was used in Secord et al. 2012. Therefore, durations given here for within-PETM bins should be considered upper estimates.

To test for potential responses within the PETM, consistent with a direct response to abiotic temperature and aridity changes, the PETM itself was subdivided into three sections: an early part of the PETM in which mean annual temperature and aridity were increasing (onset to early body), a middle part of the PETM in which mean annual temperature and aridity remained high (body), and a late part of the PETM in which climate returned to cooler and wetter background conditions (early recovery; Bowen et al. 2006; Kraus et al. 2013). To test for potential responses at the beginning and end of the PETM, specimens were binned into pre-PETM, PETM, and post-PETM groups.

Image Collection

To characterize response in terms of potential changes in tooth morphology, we first imaged specimens in 2D and 3D, so that sub-millimeter-scale features could be accurately quantified. For 3D analyses of shape and topography, all unbroken and unworn M₁s and M₂s were microcomputed tomography (microCT)-scanned on a Nikon XTH 225 ST or a GE v|tome|x M 240 (Table 2). Scan voxel resolution was similar between the two machines (3.73–7.54 μm and 5.20–6.10 μm, respectively). Data derived from CT-scanned specimens used in this study are freely available on MorphoSource (Boyer et al. 2016; Supplementary Table 1). Previously scanned, published specimens of Colpocherus and Macrocranion from the Sand Creek Divide section of the BHB were added to the dataset (Rose et al. 2012). Digital crown surfaces of each CT-scanned specimen were segmented, cropped by hand at the enamel cervical margin by a single observer (N.S.V.), and retriangulated using Avizo v. 9.1.1 (FEI Visualization Science Group, Berlin) according to the protocol suggested by a recent methodological investigation (Spradley et al. 2017). To evaluate measurement error (ME) due to cropping, the first 10% of the sample of digital crown surfaces for each of the stem erinaceid taxa were recropped a total of three times. Each of those three surfaces were also duplicated three times to measure pseudolandmark placement ME, for a total of nine copies of each specimen. For 2D analyses of both crown area and univariate features, all scanned specimens and additional specimens were photographed with a Leica EZ4 HD light microscope. Photographs were taken by a single observer (N.S.V.).

Table 2. Sample sizes of the three focal taxa for each type of measurement. M₁ sample sizes refer to specimens measured to calculate crown area. M₂ 3D sample sizes refer to specimens analyzed using 3D geometric morphometrics and dental topographic metrics. M2 1D sample sizes refer to specimens analyzed using targeted univariate measures of shape. PETM, Paleocene−Eocene thermal maximum.

Molar Size Measurement and Analyses

Molar size, and by correlation body mass, was measured by calculating the log-transformed crown area (length × width) of the M₁ (Gingerich et al. 1982; Hlusko et al. 2006; Table 2). Measurements based on photographs were collected in ImageJ by the same observer (N.S.V.). Each measurement was taken three times, and the mean of the three measurements was used in statistical analyses to reduce ME (Yezerinac et al. 1992). In addition to area, we analyzed length and width separately to see whether a potential change in crown proportions could have influenced our estimate of crown size. Centroid size, an alternative, unitless measure of size, was produced as a by-product of the 3D geometric morphometric molar shape analyses. It is usually the preferred method of quantifying size in a geometric morphometric context (Zelditch et al. 2012), but it is less strongly correlated with body mass and therefore not preferred over log-transformed crown area in this study (Gingerich and Smith 1984; St. Clair and Boyer 2016). To test hypotheses of size change between different climatic or ecological regimes, we compared time bins using a bootstrapped comparison of M₁ area, length, and width means with 1000 iterations (Kowalewski and Novack-Gottshall 2010). Based on previous morphometric studies, bins were not compared if either bin contained fewer than five specimens (Polly 2003). Pairwise tests were corrected using the false discovery rate correction (Benjamini and Hochberg 1995).

Molar Shape Measurement

The crown of M₂ was used to evaluate molar shape change based on the high heritability of its morphological features and because its central position in the molar field gives it the greatest number of tooth–tooth contacts during mastication, further constraining tooth shape for proper occlusion (Townsend et al. 2003; Table 2). The complex occlusion of the interlocking molar cusps constrains M₂ morphology to be relatively conservative due to the strong selective disadvantage of malocclusion in mammals with tall molar cusps (Gingerich and Winkler 1979).

The 3D measures chosen for this study were two dental topographic metrics and pseudolandmark-based geometric morphometrics. These were applied only to specimens with unworn, complete tooth crowns to avoid introducing statistical noise because of obliterated morphology. Three popular types of dental topographic metrics that have been shown to correlate with heuristic dietary categories in extant mammals, and are therefore one way to measure potential response in terms of tooth function, include orientation patch count and the related rotation-corrected metric (OPC, OPCR), relief index (RFI), and Dirichlet normal energy (DNE; Evans et al. 2007; Boyer 2008; Bunn et al. 2011; Pineda-Munoz et al. 2017; López-Torres et al. 2018). They are used to analyze the complete, continuous tooth crown and are not reliant on identifying points of homology. OPC and OPCR aim to approximate the number of features, or “tools” on a tooth crown (Evans et al. 2007; Pineda-Munoz et al. 2017). We did not analyze OPC or OPCR, because there was no reason to expect this kind of major morphological change in the molars of Colpocherus, Macrocranion, and Talpavoides through the studied interval. All three taxa retain the same overall morphology of six cusps lacking cingulids, cuspules, or other additional features throughout their evolutionary history (Bown and Schankler 1982; Smith et al. 2002; Beard and Dawson 2009; Rose et al. 2012). In contrast, RFI and DNE are sensitive to feature quality as well as feature quantity (King et al. 2005; Gutzwiller and Hunter 2015; Pampush et al. 2018). The RFI metric expresses crown height as a ratio of 3D occlusal area to 2D planimetric footprint area (Boyer 2008). The DNE metric integrates crown curvature across a digital surface, effectively capturing tooth sharpness (Bunn et al. 2011; Winchester et al. 2014; Pampush et al. 2016a). They were chosen because biomechanical modeling of insectivore molars predicts that cusps subject to increased risk of fracture or increased rates of wear have more robust cusps with lower relief (lower RFI) and sharpness (lower DNE). In contrast, the optimal cusps for processing invertebrates would have higher relief (higher RFI) and sharpness (higher DNE; Evans and Sanson 2003, 2005). Even without being able to predict exactly in which direction these metrics might change for a particular extinct species, we hypothesized that any number of changes in the abiotic and biotic components of the environment could have led to changes in fracture risk, rates of wear, degree of insectivory, or other functional changes that would have affected cusp relief and sharpness.

A high-throughput pseudolandmark-based geometric morphometric approach was used to rapidly characterize the shape of the entire molar crown in high detail. To measure 3D aspects of morphology, specimens were first aligned in auto3dgm in MATLAB using 256, then 2048 pseudolandmarks (Puente 2013; Boyer et al. 2015). Aligned specimens were then rotated into anatomic position using Meshlab with the occlusal plane parallel to the global x, y plane. The plyClip function in the R package molaR (Pampush et al. 2016b) was used to crop all digital crown surfaces to a common horizontal plane at the lowest point of the talonid basin, reducing ME due to cropping to <5%. RFI and DNE were then calculated using the default parameters in molaR. One specimen (USNM 538323) had to be excluded from the analysis of DNE because of surface artifacts introduced by low-resolution scanning.

The shape of the crown of M₂ was quantified using the aligned 2048-pseudolandmark configurations for each specimen output by auto3dgm. Shape error due to algorithmic pseudolandmark placement and surface cropping of each tooth crown were quantified by analyzing the replicate surfaces using previously proposed methods for evaluating error (Vitek et al. 2017). In each intraspecific dataset, pseudolandmark placement error was low (<1% ME), but cropping error had a large effect on pseudolandmark configurations even after attempted cropping standardization (ME > 15% in all datasets; Supplementary Fig. 1). Cropping-related error was clustered around the crown cervical margin, as expected, but also propagated nonrandomly through the surface, resulting in spurious amounts of variation at local minima on the crown digital model (Fig. 2A–C). To minimize the effects of this error on results, all analyses were conducted using information from only the set of 1 − N repeatable principal components (PCs) with ME of 10% or lower based on a principal components analysis (PCA) of the dataset including replicate specimens (repeatable PCs: Colpocherus N = 3, Macrocranion N = 4, Talpavoides N = 3; Yezerinac et al. 1992; Fruciano 2016; Vitek et al. 2017).

Figure 2. M₂ crown shapes. A−C, Maps of where shape variation due to cropping error are concentrated on the crown surface, in replicates of single specimens as examples: A, UF 328055, Colpocherus; B, UF 326899, Macrocranion; C, UF 327009, Talpavoides. Note that this spurious variation is not limited to the cervical margin of the crown, where the error was introduced. D−I, Maps of each species showing patterns of differences between samples of specimens. Crown shapes predicted by the minima (D, F, H) and maxima (E, G, I) of PC 1 for Colpocherus (D, E), Macrocranion (F, G), and (H, I) Talpavoides (H, I). J, K, Mean differences between crown shapes of the early and mid-Paleocene–Eocene thermal maximum (PETM) mapped onto mean early-PETM shape in Colpocherus (J) and Macrocranion (K). Warmer colors indicate either greater differences (D−K) or greater variability (A−C) within each 2048-pseudolandmark dataset analyzed. Areas of cool colors and mottled patterns in relatively cool colors can indicate spurious or nonsignificant variation. met, metaconid; taln, talonid notch; triba, trigonid basin. (Color online.)

After error was evaluated and replicate surfaces were removed, 3D shape was explored, starting with a PCA to summarize the primary structure of shape variation within each species. To visualize the shape differences between pairs of successive stratigraphic bins, heat maps were constructed using Meshlab and custom R code. Visual comparison of the heat maps overlaid on the pair of mean shapes was used to develop hypotheses of any possible changes that could have occurred between two bins (e.g., Fig. 3). We then developed a suite of univariate shape descriptors calculated from linear and angular measurements of crown features to test each of those hypotheses (Fig. 4). Those descriptors could be used to quantitatively test for differences in features that appeared qualitatively different in 3D heat maps, because the descriptors could be measured in a larger sample of broken and worn specimens that were not suitable for 3D analyses of the complete crown (Table 2). Each linear and angular measurement was done three times by a single observer (N.S.V.), and the mean of each trio of measurements was used in downstream analyses to minimize ME (Yezerinac et al. 1992).

Figure 3. Shape differences between time bins in the M₂ of Macrocranion (A−D) and Talpavoides (E–L). Between-bin shape for Macrocranion: mean Paleocene–Eocene thermal maximum (PETM) shape (A, C) versus mean post-PETM shape (B, D) in occlusal (A, B) and posterior (C, D) views. Between-bin shape for Talpavoides: mean PETM shape (E, G) versus mean pre-PETM shape (F, H) in occlusal (E, F) and posterior (G, H) views; mean PETM shape (I, K) versus mean post-PETM shape (J, L) in occlusal (I, J) and posterior (K, L) views. All shapes colored by the differences between the two bins compared in a row, with hotter colors corresponding to greater difference, as in Fig. 2. co, cristid obliqua; hyp, hypoconid; hypc, hypoconulid; proco, protoconid; triba, trigonid basin; talba, talonid basin. (Color online.)

Figure 4. Illustration of the univariate measurements taken post hoc on M₂s. An M₂ of Macrocranion, UF 283308, is shown as an example in occlusal (A), lingual (B), and posterior (C) views. CA, crown area; HHID, hypoconid−hypoconulid intercusp distance; L, length; MEID, metaconid−entoconid intercusp distance; ML, metaconid length; R-, relative; TH, trigonid height; TW, talonid width; W, width.

We did not hypothesize any functional meaning or particular cause for these descriptors. Their only source was the heat maps, and their purpose was to complement analyses of 3D shape that were more strongly limited by sample size than analyses of univariate shape metrics. Only after this analysis was completed did we evaluate whether any observed, statistically significant changes in univariate measures of shape could be explained by function related to diet. To do so, we investigated whether changes occurred on functional surfaces that are arguably subject to selection related to mechanical requirements for chewing food with different structural qualities (Ungar et al. 2018). The locations of changes on the tooth crown were compared with the order of appearance and distribution of wear facet development to determine whether differences in dental topography are related to changes in chewing function (Kay and Hiiemae 1974). Wear facets form through chewing, and if the locations of changes were also the locations of wear facets, we inferred that the morphological changes had some relationship to functional changes during chewing (Green and Croft 2018). We used order of appearance to infer which facets formed more quickly than others.

Molar Shape Analyses

To statistically test for differences in overall 3D crown morphology in different climatic or ecological regimes, differences in 3D pseudolandmark configurations between bins were evaluated using Procrustes analysis of variance (ANOVA) in the geomorph package (Adams et al. 2017). A minimum sample size of N = 5 was imposed to test a bin for significant differences (Polly 2003). For dental topographic measures and univariate measures of shape, the presence of significant differences between bins was tested using a bootstrapped comparison of means with 1000 iterations (Kowalewski and Novack-Gottshall 2010). Pairwise tests were corrected using the false discovery rate correction (Benjamini and Hochberg 1995).

All analyses were performed in R v. 3.5.3 (R Core Team 2015). Live versions of the R code used to conduct analyses and produce visualizations are stored on GitHub under the project name “PETM-erinaceomorphs” (https://github.com/nsvitek/PETM-erinaceomorphs). An archival snapshot of the code, supplemental tables of analytical results, the pseudolandmark configurations output by auto3dgm and a table of specimen numbers, MorphoSource media groups, and metadata are reposited on the Dryad Digital Repository (https://doi.org/10.5061/dryad.mkkwh70zr). A list of voucher specimens used in this study is provided in Supplementary Table 2.

Results

First and Last Appearances within the BHB

Before this survey, biostratigraphic extents of these genera were incompletely known for this interval, with Talpavoides unreported for Wa-0, the temporal range of Colpocherus unknown outside of two localities from the body of the PETM in the Sand Creek Divide area of the BHB, and replacement of Macrocranion junnei by Macrocranion nitens proposed to occur at the Wa-0/1 boundary (Smith et al. 2002; Yans et al. 2006; Rose et al. 2012). Based on our survey of 285 total cataloged specimens from the pre-PETM, 12,664 cataloged specimens from the PETM, and 2454 cataloged specimens from the post-PETM, Colpocherus is only present within the PETM during Wa-0, Macrocranion junnei (hereafter, Macrocranion) is present from the PETM to post-PETM during Wa-0–Wa-2, and Talpavoides is present throughout the studied interval from pre- to post-PETM during Cf-3–Wa-2.

Size Change Limited to PETM Boundaries

The M₁ crown area of Colpocherus did not change significantly from its first appearance in the early PETM to the middle PETM (N = 9, 20; p = 0.75; Fig. 1C). Similarly, the M₁ crown area of Macrocranion did not change from the early to mid-PETM (N = 36, 19; p = 0.52; Fig. 1D), but it did significantly increase between the PETM and post-PETM (N = 56, 7; p = 0.002; Fig. 1D). Crown area of Talpavoides did not clearly change within or between any time bin evaluated here (N = 6, 15, 27; p > 0.21; Fig. 1E, Supplementary Tables 3,4). Length and width measurements of the M₁ both showed the same patterns as area measurements (Supplementary Fig. 2, Supplementary Tables 3,4).

Shape Change Limited to PETM Boundaries

The PCA of all 2048 pseudolandmarks on M₂ shows overlap in morphology between specimens of different time bins. This is true of all taxa studied (Supplementary Fig. 3). PC 1 described some differences around the crown cervical margin and local minima of the occlusal surface (Fig. 2D–I, buccal view, lingual view) as well as regions on the buccal half of the tooth in Macrocranion and Talpavoides (Fig. 2D–I, occlusal view). Shape differences between within-PETM bins as measured by Procrustes ANOVA of repeatable PCs were not statistically clear for any taxon (5 < N < 10; p > 0.15; Supplementary Table 3). Differences in mean shape between subdivisions of the PETM were most prominent around the trigonid basin and talonid notch (Fig. 2J,K).

Although PETM and post-PETM specimens of Macrocranion do not occupy different regions of 3D pseudolandmark morphospace (p = 0.955; Supplementary Table 4), the mean shape differences between PETM and post-PETM bins were not randomly distributed across the tooth (Fig. 3A–D). The differences that did not correspond to patterns of shape variation due to cropping error were consistent with a change in molar canting angle, or the degree to which the cusps tilt lingually relative to the base of the tooth. These differences were not apparent to the authors before the 3D geometric morphometric analyses and motivated further investigation.

When mean shapes of the pre- and post-PETM bins of Talpavoides were compared with that of the PETM bin, differences were concentrated on the protoconid, the cristid obliqua, and between the hypoconid and hypoconulid (Fig. 3E–L). Visualized differences between pre-PETM and PETM bins additionally corresponded to a change in buccolingual curvature of the cusps. Similar to Macrocranion, these described differences were not obvious before 3D geometric morphometric analyses.

To better test the biological reality of the differences between the average shapes of bins with relatively low 3D sample sizes, we developed and collected univariate measurements that targeted proposed aspects of shape differences. In Colpocherus, two characters, the length of the metaconid relative to the total M₂ length and the distance between the metaconid and entoconid, are predicted by higher levels of change between early- and mid-PETM surfaces. Potential changes in both features lack statistical clarity when measured directly on a larger sample size (p > 0.196; Dushoff et al. 2019; Fig. 5A, Supplementary Table 3). Three-dimensional dental topographic measures also lacked significant change between bins (N = 7, 6; p > 0.171; Supplementary Table 2, Supplementary Fig. 3).

Figure 5. Univariate measurements developed from qualitative between-bin differences in Figs. 3 and 4 compared among sampling bins for each species: A, Colpocherus, B, Macrocranion, and C, Talpavoides. Dashed lines indicate the beginning and end of the carbon isotope excursion delimiting the Paleocene–Eocene thermal maximum (PETM). Plotting shapes indicate biozones. Colors indicate sampling bins used in this study. RHHID, relative hypoconulid−hypoconid intercusp distance. Light brackets indicate bins being compared. Heavier brackets indicate statistical comparisons. Asterisks next to brackets indicate statistically significant differences. (Color online.)

In Macrocranion, no difference between early- and mid-PETM surfaces could be found that did not also correspond to cropping error, and no additional measurements were collected. Canting angle decreased significantly between the PETM and post-PETM, as predicted by differences between 3D shapes (N = 75, 10; p = 0.001; Figs. 4A–D, and 5B, Supplementary Table 4), but not within the PETM, where no such change was predicted (N = 35, 37; p = 0.231; Supplementary Table 3). A pattern of minimal change was observed in diet-related dental topographic metrics RFI and DNE that could not be tested for significance due to sample size constraints (Supplementary Fig. 4). The surfaces most affected by decreased canting angle correspond with those that are more quickly affected by tooth wear during the chewing cycle compared with other surfaces (Supplementary Fig. 5).

In Talpavoides, changes across PETM boundaries appear concentrated in the trigonid basin, as measured by relative trigonid height, and in the position of the entoconid, as evaluated by both relative talonid width and the relative distance between the hypoconid and hypoconulid (RHHID; Fig. 4). The relative hypoconid–hypoconulid distance increased significantly from the PETM to the post-PETM (N = 12, 26; p = 0.002; Fig. 5C, Supplementary Table 4), but not within the PETM (N = 5, 7; p = 0.542; Supplementary Table 3), where changes could not be predicted due to a lack of complete specimens from the mid-PETM. This measurement corresponds to the length of the postcristid and affects the space available to develop facet 4 of Kay and Hiiemae (1974), one of the first wear facets to form on the teeth of Talpavoides (Supplementary Fig. 5). After correcting for multiple tests, the relative height of the trigonid and relative width of the talonid were found to be statistically unclear (Supplementary Tables 3,4). The statistically clear shape changes in Talpavoides do not appear to correspond to changes in the 3D topographic metrics RFI and DNE. Note that this qualitative observation could not be quantitatively tested for significance due to sample size constraints (Supplementary Fig. 4).

Discussion

Timing and Function of Change

The timing of morphological change within a species is the primary way to discriminate abiotically from biotically driven changes in the predictive framework that we developed for the PETM. Specifically, morphological changes within the PETM, when abiotic climate continued to change but the community remained relatively stable, are only consistent with an explanation of abiotic drivers. Patterns of static morphology within the PETM are more consistent with the primary influence of biotic drivers. Both biotically and abiotically driven changes could occur at the boundaries of the PETM, when both abiotic climate and biotic community composition were changing.

Our results do not support a hypothesis of abiotic change as a direct driver of altered dental morphology in these three species of erinacemorphs. The crown size of M₁ and crown surface topography of M₂ appear to be static within the PETM in all three taxa, even during the extensive climatic change documented in the same or correlated sections in the BHB (Clyde and Gingerich 1998; Gingerich 2006; McInerney and Wing 2011; Wing and Currano 2013). Although 3D metrics did not help discriminate between the two drivers, it is possible that the lack of significant change in 3D metrics is due to a low sample size in some cases. In particular, Talpavoides was recovered too rarely within the PETM to evaluate for significant change (Table 2), and we cannot evaluate consistency with either hypothesis for this species. A more complete suite of complementary 3D and 2D results for Macrocranion and Colpocherus is consistent with a true lack of change during the PETM. To evaluate whether lack of statistical significance in 3D metrics masked true morphological change, heat maps were constructed to visualize the difference between mean early- and mid-PETM molar shapes for each species. The potential changes to tooth morphology suggested by patches of greater change on those heat maps were largely consistent with differences related to cropping error, particularly local minima on the crown surfaces such as the talonid notch and trigonid basin. The potential for real biological differences on other parts of the tooth was not upheld by targeted measures of morphology taken on a much larger sample of specimens (Table 2), consistent with the idea that those potential changes are more likely the spurious result of cropping error or small sample sizes, not real biological change. All of the changes predicted by pseudolandmark differences between early- and mid-PETM shapes lacked statistical clarity in Colpocherus.

Because the hypothesis of biotically driven response is the same as our null hypothesis in this study, it was not tested in the same way as the alternative explanation, and we cannot say that our results definitively support this hypothesis (Hunt et al. 2008). However, overall, the timing of changes observed is consistent with explanations of biotically driven, community-mediated responses to climate changes. This stands in contrast with prior investigations of PETM fauna that have suggested a primary role for abiotic drivers in producing morphological change (e.g., equid body size: Secord et al. 2012). In the future, likelihood models may allow for a more balanced statistical comparison of the abiotic and biotic response hypotheses.

The morphological patterns of stem erinaceid species that cross the end-PETM boundary in the BHB are nonetheless interesting in their own right, even though they do not help distinguish between biotic and abiotic drivers. If we had observed no morphological change across the PETM boundaries, these species would have been a remarkable example of stasis and potential resilience in the face of extreme climate change (Bowen et al. 2006). As it is, most aspects of dental morphology appear to remain static across the boundary in both of these species, Macrocranion and Talpavoides. Size changes only in Macrocranion (Fig. 1). Shape remained unchanged across much of the crown surface, as indicated by large swaths of cool colors interrupted by localized regions of change on heat maps (Fig. 5). Low sample size in one or more bins precluded evaluating whether these localized changes would have resulted in clear, significant changes to any 3D metric. However, 3D analyses still contributed useful information, in that the heat maps identified those localized changes that we had not previously considered. Not all of those potential changes were upheld by results from larger samples of 2D measurements (Supplementary Table 4), consistent with the overall pattern of only small amounts of change concentrated in only a few features of the tooth.

In both species that cross the end-PETM boundary, those small components of shape change are localized on the buccal side of the molar, particularly on the cusps and crests that are involved in shearing functions during the chewing cycle (Fig. 3). These morphological changes, though only affecting small subsections of the molar crown surface, were likely functional, and from that perspective were likely under selection. If environmental changes selected on aspects of molar morphology through shifting diets, then changes should be observable in features correlated with comminution of food with different material properties or on functional surfaces used in food processing (Kay and Hiiemae 1974; Kay 1975; Janis and Fortelius 1988; Evans et al. 2007; Boyer 2008). Although we lacked a sufficient number of complete tooth crowns to test for end-PETM changes in the dental topographic metrics DNE and RFI, we could evaluate morphological change with regard to wear facets. Wear facets are a direct record of which parts of the occlusal surface are interacting during the chewing cycle (Crompton 1971; Kay and Hiiemae 1974). Greater amounts of wear or earlier appearance of wear facets implies that those regions of the tooth see greater use than others (Fortelius 1985).

Molar shape changes in both Macrocranion and Talpavoides occur primarily on buccal crests and cusps (Fig. 3). Early-developing wear facets on the lower molars of those same species also form on the buccal crests and cusps. In Macrocranion, changes in canting angle affect some of first occlusal surfaces to be used and worn down in molar function (Fig. 5, Supplementary Fig. 5). The wear facets are consistent with shearing functions during phase I of chewing, producing facets 1, 3, and 4 within the framework developed to describe primate mastication and tooth wear by Kay and Hiiemae (1974) and with shearing facets 2, 3, and 6 of the marsupial Didelphis marsupialis in the framework of Crompton and Hiiemae (1970). The reduction in canting angle increases the height of the metaconid, entoconid, and cristid obliqua relative to the bottom of the talonid basin, especially those parts available to act as shearing surfaces. By increasing the height of cusp and crest available for function, the change in canting angle may have given the teeth of Macrocranion a longer functional life span or may indicate an increase in the amount of tough, ductile food being ingested (Kay and Hiiemae 1974; Strait 1993). In Talpavoides, changes in intercusp distances correspond to a shortening of the postcristid after the PETM. Evidence of early wear facet formation (Supplementary Fig. 5) along the postcristid suggests it was in heavy use during Talpavoides mastication. This wear facet, along with those that develop on the protocristid (Facet 1) and cristid obliqua (Facet 3), is consistent with shearing functions during chewing (Kay and Hiiemae 1974). The postcristid functions as a shearing surface, and its reduced size is consistent with reduced toughness and ductility and increased brittleness of foods eaten by small, insectivorous mammals, or a relatively higher proportion of “hard-bodied” insects such as beetles (Strait 1993). This correspondence supports the hypothesis that shape changes are localized on functionally important features and are therefore likely targets of selection (Barrett and Hoekstra 2011). However, potential functional interpretations of these facets require further development, because changes in mesowear such as these can also be due to changes in environmental factors such as the amount of exogenous grit as well as changes in diet (Loffredo and DeSantis 2014; Green and Croft 2018).

Study Limitations

One potential bias to our conclusions is the possibility that we did not fully consider all aspects of morphological change in the molars under study and might have failed to observe changes that occur within the PETM. However, certain exploratory steps were included in our analytical protocol specifically to distinguish the absence of evidence for change from the evidence for absence of change. In particular, we calculated heat maps of differences between average tooth shapes of bins to exhaustively search for potential changes, then tested those potential changes for statistical significance, and by proxy biological reality, using targeted univariate metrics that captured the potential morphological change. Based on the results that the heat maps included potential morphological changes that were not upheld by further measurement, we conclude that the heat-map approach behaved relatively liberally in its characterization of potential shape change, was not an overly conservative exploratory tool, and adequately captured all potential change on the M₂ crown. Therefore, we argue that we comprehensively characterized dental shape as well as is reasonably possible with the tools currently available.

Other potential limitations of our study apply widely to almost all studies of the fossil record. Based on our results, we cannot definitively say that these species did not respond at all to abiotic climate change in the PETM. It is possible that other anatomic systems, including shape of the M₁, size of the M₂, morphology of the M₃, or additional teeth and bones, could contain other patterns of change through time. Traits more difficult to document in the fossil record, such as fur density, behavior, or life history, may have been much more sensitive to temperature change than molar morphology (Storz et al. 2019). However, previous studies of molar morphology have found that it can closely track isotopically inferred changes in diet (Kimura et al. 2013). Compared with other traits measured using similar time intervals to those used in this study (Table 1, compare with 1.45–84.5 kyr intervals of Clyde and Gingerich [1994] and Wood et al. [2007]), molar size and shape have similar rates of evolution, so we do not consider molar morphology to be a set of especially insensitive traits (Clyde and Gingerich 1994; Wood et al. 2007; Uyeda et al. 2011). Study of additional systems and the potential for mosaic evolution is an interesting avenue of future research, but one that will require considerably more fossil collection, measurement, and analysis for the focal taxa in this study. Other potential limitations of the study are more difficult to account for, as for any study of the fossil record. For example, it is possible that an unmeasured abiotic variable, such as maximum summer temperature or minimum winter temperature, changed in a pattern much more similar to the biotic variables considered, unlike the abiotic variables we were able to account for. In addition, we cannot discriminate between specific abiotic or specific biotic drivers in this study. Sample size also acts as a current limit on many potentially interesting avenues of study, despite our efforts to choose the most common species from a large collection of fossils gathered through intensive sampling efforts over many field seasons. Even with all of these limitations in mind, the consistency of observed patterns with biotically driven change, especially within the PETM, when many of these limitations are minimized, is a remarkable contrast to contemporaneous change in other species, such as equids (Secord et al. 2012). For tooth shape within these species not to change within the PETM although it did at other times and for tooth shape not to change within the PETM as climate continued to shift highlight that responses to climate change may not be driven by a linear relationship to abiotic climate alone.

Summary and Future Directions

Overall, the three stem erinaceid species studied here appear to be largely insensitive to climate change, even when it is related to clear shifts in biome, in both dental morphology and reconstructed body size. That pattern differs from the patterns in other previously documented taxa that seem more influenced by mean annual temperature, such as the Pleistocene woodrat Neotoma (Smith et al. 1995) and equid body size during the PETM (Secord et al. 2012). It also raises questions of what proximate agents of selection on each documented insectivore phenotype remained relatively unchanged during the PETM. Different ecological functional groups are primarily limited in realized niche by different factors (Holt 1990). For some, like the herbivorous Pleistocene Neotoma and PETM equids, physiological temperature limitations may be the primary limitation on body size. For those species, climate changes may more directly impact the selective pressures on their populations. For others, competition for the same limited resource among a large number of species may limit adaptation to other environmental changes and be the primary control on response (Johansson 2008). Given the strong ecological structuring and mixed evidence for physiological ties between body size and temperature in extant insectivorous mammal communities, lack of response to abiotic climate change in extinct insectivores is not necessarily unexpected (Churchfield et al. 1999; Ochocinska and Taylor 2003; Yom-Tov and Yom-Tov 2005; White and Searle 2007; but see Contoli et al. [2000] and Rychlik et al. [2006] for contrary expectations). Additional study of medium-sized herbivorous mammals, such as the “condylarths” Ectocion and Copecion (Gingerich 2006), as well as other small-sized insectivores during the PETM, should be performed to test the hypothesis that larger herbivorous species are more likely to respond directly to abiotic change, while small-bodied species are more likely to be primarily limited by their community ecological context.

As a test case, results of this study are not yet widely applicable to other taxa. Instead, the results demonstrate proof of concept for the research framework and an intriguing contrast to the responses of other mammals to climate change (Smith et al. 1995; Secord et al. 2012). Based on the range of responses thus far documented in the PETM, we suggest that body size may help to predict whether abiotic climate change is likely to produce direct or community-mediated effects on a species. Very small species (body mass <50 g) may be more affected by biotic change, on average, than larger species. Smaller-bodied mammalian species may also be less noticeably affected than larger ones, similar to patterns documented in response to current anthropogenic climate change (Purvis et al. 2000; Cardillo and Bromham 2001; McCain and King 2014). For example, of the three very small-bodied genera studied here and the eight additional genera of very small mammals for which body size has been estimated in the PETM of the BHB, only two have smaller PETM representatives than pre- or post-PETM representatives (Secord et al. 2012; Felibert et al. 2017). In contrast to the 40% of the rest of the mammalian fauna that shows a PETM-related shift toward smaller body mass, only 18% of the very small-bodied genera documented so far have a similar response. Variation in which factors ultimately limit body size may also explain why some small mammals are less likely to conform to patterns more apparent in larger mammals, such as Bergmann's rule (Meiri and Dayan 2003; Smith et al. 2004). It is possible that the rules of response are different for small-bodied mammals, which have different energy constraints and physiological tolerances in comparison to large-bodied mammals (Schmidt-Nielsen 1984; Terry and Rowe 2015). A better understanding of how ecological variables influence adherence to biological rules and relationships to climate may improve our ability to predict mammalian responses to modern climate change (Purvis et al. 2000).