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. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013). 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. Reference Gill, Williams, Jackson, Lininger and Robinson2009; van Gils et al. Reference Gils, Lisovski, Lok, Meissner, Ożarowska, de Fouw, Rakhimberdiev, Soloviev, Piersma and Klaassen2016). 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. Reference Gill, Williams, Jackson, Lininger and Robinson2009). Examples of such fossil records are rare (Gill et al. Reference Gill, Williams, Jackson, Lininger and Robinson2009; Terry Reference Terry2018). 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 Reference Strait2001; Wing et al. Reference Wing, Harrington, Smith, Bloch, Boyer and Freeman2005; Yans et al. Reference Yans, Strait, Smith, Dupuis, Steurbaut and Gingerich2006; Kraus and Riggins Reference Kraus and Riggins2007; Smith et al. Reference Smith, Wing and Freeman2007; Chester et al. Reference Chester, Bloch, Secord and Boyer2010; Adams et al. Reference Adams, Kraus and Wing2011; Rose et al. Reference Rose, Chester, Dunn, Boyer and Bloch2011, Reference Rose, Chew, Dunn, Kraus, Fricke and Zack2012; Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012; Baczynski et al. Reference Baczynski, McInerney, Wing, Kraus, Bloch, Boyer, Secord, Morse and Fricke2013, Reference Baczynski, McInerney, Wing, Kraus, Bloch and Secord2017; Kraus et al. Reference Kraus, McInerney, Wing, Secord, Baczynski and Bloch2013; Bourque et al. Reference Bourque, Howard Hutchison, Holroyd and Bloch2015; Morse et al. Reference Morse, Chester, Boyer, Smith, Smith, Gigase and Bloch2019). 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. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012; Kraus et al. Reference Kraus, McInerney, Wing, Secord, Baczynski and Bloch2013; Baczynski et al. Reference Baczynski, McInerney, Freeman and Wing2019). 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. Reference Koch, Zachos and Gingerich1992; Gingerich Reference Gingerich2001; Wing et al. Reference Wing, Harrington, Smith, Bloch, Boyer and Freeman2005). 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 Reference Kurtén1960; Smith et al. Reference Smith, Betancourt and Brown1995; McGuire Reference McGuire2010; Prothero et al. Reference Prothero, Syverson, Raymond, Madan, Molina, Fragomeni, DeSantis, Sutyagina and Gage2012; Williams and Blois Reference Williams and Blois2018). However, climate change affects many aspects of an ecosystem simultaneously (Davis et al. Reference Davis, Jenkinson, Lawton, Shorrocks and Wood1998; Alexander et al. Reference Alexander, Diez and Levine2015; Eskelinen et al. Reference Eskelinen, Kaarlejärvi and Olofsson2017; Mokany et al. Reference Mokany, Bush and Ferrier2019). 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. Reference Davis, Jenkinson, Lawton, Shorrocks and Wood1998; Cahill et al. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013). This range of potential proximate drivers introduces complexity into predictions of how any single species will ultimately respond to climate change (Graham Reference Graham, Diamond and Case1986; Roy and Pandolfi Reference Roy, Pandolfi, Lovejoy and Hannah2005; Stewart Reference Stewart2008). 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. Reference Gill, Williams, Jackson, Lininger and Robinson2009; Cahill et al. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013; van Gils et al. Reference Gils, Lisovski, Lok, Meissner, Ożarowska, de Fouw, Rakhimberdiev, Soloviev, Piersma and Klaassen2016).
Proximate drivers can be broadly categorized as abiotic and biotic (Cahill et al. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013). 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. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013). 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 Reference Smith and Betancourt2006), 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. Reference Smith, Betancourt and Brown1995).
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. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013). 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 Reference Terry2018). Biotic factors are a frequent cause of local population declines that are ultimately related to climate change (Cahill et al. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013; Westover and Smith Reference Westover and Smith2020).
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. Reference Alexander, Diez and Levine2015; Eskelinen et al. Reference Eskelinen, Kaarlejärvi and Olofsson2017). Both categories of drivers can ultimately be driven by change in abiotic climate (Cahill et al. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013). 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 Reference Gingerich2006; McInerney and Wing Reference McInerney and Wing2011). The PETM occurred at the beginning of the Eocene epoch, ~56 Ma (Gradstein et al. Reference Gradstein, Ogg and Hilgen2012). 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 Reference Clyde and Gingerich1998; Gingerich Reference Gingerich2006; Woodburne et al. Reference Woodburne, Gunnell and Stucky2009; McInerney and Wing Reference McInerney and Wing2011; Wing and Currano Reference Wing and Currano2013).
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. Reference Bowen, Bralower, Delaney, Dickens, Kelly, Koch, Kump, Meng, Sloan, Thomas, Wing and Zachos2006; Aziz et al. Reference Aziz, Hilgen, van Luijk, Sluijs, Kraus, Pares and Gingerich2008; Westerhold et al. Reference Westerhold, Röhl, McCarren and Zachos2009; McInerney and Wing Reference McInerney and Wing2011). 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. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012; Fig. 1A).
Patterns of M1 crown area through the Paleocene–Eocene thermal maximum (PETM) (A) in comparison to the oxygen isotope record reconstructing mean annual temperature from δ18O of tooth enamel in the mammal Coryphodon. B, Patterns of M1 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 δ18O values. A, B, Data from Secord et al. (Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012). 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. Reference Foreman, Heller and Clementz2012; Kraus et al. Reference Kraus, McInerney, Wing, Secord, Baczynski and Bloch2013; Baczynski et al. Reference Baczynski, McInerney, Wing, Kraus, Bloch and Secord2017). The timing of paleohydrologic changes largely follows the timing of paleotemperature changes (Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012; Kraus et al. Reference Kraus, McInerney, Wing, Secord, Baczynski and Bloch2013; Baczynski et al. Reference Baczynski, McInerney, Freeman and Wing2019). During the PETM onset and into the early PETM body, climate in the BHB became warmer and more arid (Kraus et al. Reference Kraus, McInerney, Wing, Secord, Baczynski and Bloch2013). It remained arid with more seasonal or episodic rainfall through the mid-body of the PETM (Wing et al. Reference Wing, Harrington, Smith, Bloch, Boyer and Freeman2005; Foreman et al. Reference Foreman, Heller and Clementz2012; Kraus et al. Reference Kraus, McInerney, Wing, Secord, Baczynski and Bloch2013; Baczynski et al. Reference Baczynski, McInerney, Wing, Kraus, Bloch and Secord2017, Reference Baczynski, McInerney, Freeman and Wing2019). 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 Reference Barnosky2001). 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. Reference Wing, Harrington, Smith, Bloch, Boyer and Freeman2005; Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012; Baczynski et al. Reference Baczynski, McInerney, Wing, Kraus, Bloch, Boyer, Secord, Morse and Fricke2013, Reference Baczynski, McInerney, Wing, Kraus, Bloch and Secord2017; Kraus et al. Reference Kraus, McInerney, Wing, Secord, Baczynski and Bloch2013). A high-resolution stratigraphic framework (Baczynski et al. Reference Baczynski, McInerney, Wing, Kraus, Bloch, Boyer, Secord, Morse and Fricke2013) 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. Reference Baczynski, McInerney, Wing, Kraus, Bloch, Boyer, Secord, Morse and Fricke2013). In several cases, these efforts produced sample sizes that are large enough to quantify simultaneous responses in multiple taxa (Polly Reference Polly2003; Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012). 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. Reference Reed, Kaufman and Kaufman2006; Eronen et al. Reference Eronen, Evans, Fortelius and Jernvall2010; Terry et al. Reference Terry, Li and Hadly2011). 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. Reference Baczynski, McInerney, Wing, Kraus, Bloch, Boyer, Secord, Morse and Fricke2013).
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. Reference Holroyd, Hutchison and Strait2001; Gingerich Reference Gingerich2006; Smith Reference Smith2009; Smith et al. Reference Smith, Hasiotis, Kraus and Woody2009; Wing and Currano Reference Wing and Currano2013; Bourque et al. Reference Bourque, Howard Hutchison, Holroyd and Bloch2015). 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 Reference Gingerich1989; Strait Reference Strait2001; Smith et al. Reference Smith, Bloch, Strait and Gingerich2002; Heinrich et al. Reference Heinrich, Strait and Houde2008; Chester et al. Reference Chester, Bloch, Secord and Boyer2010; Rose et al. Reference Rose, Chew, Dunn, Kraus, Fricke and Zack2012). Biostratigraphic and biogeographic studies have also been undertaken (e.g., Gingerich Reference Gingerich2001; Bowen et al. Reference Bowen, Clyde, Koch, Ting, Alroy, Tsubamoto, Wang and Wang2002; Gingerich and Smith Reference Gingerich and Smith2006; Smith et al. Reference Smith, Rose and Gingerich2006; Chew Reference Chew2009; Solé and Smith Reference Solé and Smith2013; Rankin et al. Reference Rankin, Fox, Barrón-Ortiz, Chew, Holroyd, Ludtke, Yang and Theodor2015; Morse et al. Reference Morse, Chester, Boyer, Smith, Smith, Gigase and Bloch2019; Fraser and Lyons Reference Fraser and Lyons2020; van der Meulen et al. Reference van der Meulen, Gingerich, Lourens, Meijer, van Broekhuizen, van Ginneken and Abels2020), 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. Reference Koch, Zachos and Gingerich1992; Gingerich Reference Gingerich2001; Wing et al. Reference Wing, Harrington, Smith, Bloch, Boyer and Freeman2005; although see Hooker Reference Hooker2015). 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. Reference Currano, Laker, Flynn, Fogt, Stradtman and Wing2016). The plant and insect communities turned over again through the PETM recovery (Wing and Currano Reference Wing and Currano2013). By the end of the PETM, incidence of insect damage was similar to that of pre-PETM levels (Currano et al. Reference Currano, Laker, Flynn, Fogt, Stradtman and Wing2016). 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. Reference Bowen, Koch, Gingerich, Bains and Corfield2001; Gingerich Reference Gingerich2001; Gingerich and Smith Reference Gingerich and Smith2006). 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 Reference Chew2009). 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. Reference Lavergne, Mouquet, Thuiller and Ronce2010; Walther Reference Walther2010; Willis et al. Reference Willis, Bailey, Bhagwat and Birks2010), 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 Reference Clyde and Gingerich1994; Wood et al. Reference Wood, Zelditch, Rountrey, Eiting, Sheets and Gingerich2007). 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. Reference Foreman, Heller and Clementz2012; Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012; Kraus et al. Reference Kraus, McInerney, Wing, Secord, Baczynski and Bloch2013; 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 Reference Soberón and Peterson2005) and the use of vertebrate morphology as paleoclimate proxies (Dayan et al. Reference Dayan, Simberloff, Tchernov and Yom-Tov1991; Head et al. Reference Head, Bloch, Hastings, Bourque, Cadena, Herrera, Polly and Jaramillo2009; McGuire Reference McGuire2010), as well as results from previous research on the topic (Smith et al. Reference Smith, Betancourt and Brown1995; Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012), predict that such a direct relationship between abiotic climate change and species morphological responses exists (Polly et al. Reference Polly, Eronen, Fred, Dietl, Mosbrugger, Scheidegger, Frank, Damuth, Stenseth and Fortelius2011; Lawing et al. Reference Lawing, Eronen, Blois, Graham and Polly2017).
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. Reference Barnosky, Bell, Emslie, Goodwin, Mead, Repenning, Scott and Shabel2004; Terry and Rowe Reference Terry and Rowe2015). 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 Reference Yom-Tov and Yom-Tov2005; White and Searle Reference White and Searle2007). 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. Reference Novacek, Bown and Schankler1985; Gunnell et al. Reference Gunnell, Bown, Hutchison, Bloch, Janis, Gunnell and Uhen2008; Beard and Dawson Reference Beard and Dawson2009; He et al. Reference He, Chen, Gould, Yamaguchi, Ai, Wang, Zhang and Jiang2012). Postcranial fossils from the BHB support the interpretation that the three studied species were terrestrial or scansorial (Penkrot et al. Reference Penkrot, Zack, Rose, Bloch, Sargis and Dagosto2008; Manz et al. Reference Manz, Chester, Bloch, Silcox and Sargis2015; Penkrot and Zack Reference Penkrot and Zack2016). Eulipotyphlan morphology can track changes in temperature, precipitation, and food availability on a decadal scale (Yom-Tov and Yom-Tov Reference Yom-Tov and Yom-Tov2005; Poroshin et al. Reference Poroshin, Polly and Wójcik2010). 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 Reference Kay1975; Gingerich et al. Reference Gingerich, Smith and Rosenberg1982; Legendre Reference Legendre1986; Hlusko et al. Reference Hlusko, Lease and Mahaney2006; Pineda-Munoz et al. Reference Pineda-Munoz, Lazagabaster, Alroy and Evans2017; Thiery et al. Reference Thiery, Guy and Lazzari2017) that are fundamental aspects of species ecology and are likely to be perturbed by climate change (Bergmann Reference Bergmann1848; Schmidt-Nielsen Reference Schmidt-Nielsen1984; Cahill et al. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013). Body size fundamentally constrains physiology and both biotic and abiotic interactions (Calder Reference Calder1983; Schmidt-Nielsen Reference Schmidt-Nielsen1984). Rapid shifts in body size that correspond with intervals of climate change have been empirically documented (Smith et al. Reference Smith, Betancourt and Brown1995; Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012; Teplitsky and Millien Reference Teplitsky and Millien2014; van Gils et al. Reference Gils, Lisovski, Lok, Meissner, Ożarowska, de Fouw, Rakhimberdiev, Soloviev, Piersma and Klaassen2016). Dietary resources are the primary way that organisms interact with the environment in terms of their dental morphology (Kay Reference Kay1975; Sheine and Kay Reference Sheine and Kay1977; Gingerich et al. Reference Gingerich, Smith and Rosenberg1982; Coiner-Collier et al. Reference Coiner-Collier, Scott, Chalk-Wilayto, Cheyne, Constantino, Dominy, Elgart, Glowacka, Loyola, Ossi-Lupo, Raguet-Schofield, Talebi, Sala, Sieradzy, Taylor, Vinyard, Wright, Yamashita, Lucas and Vogel2016). Food availability is an important biotic cause of species responses to climate change (Cahill et al. Reference Cahill, Aiello-Lammens, Fisher-Reid, Hua, Karanewsky, Yeong Ryu, Sbeglia, Spagnolo, Waldron, Warsi and Wiens2013; van Gils et al. Reference Gils, Lisovski, Lok, Meissner, Ożarowska, de Fouw, Rakhimberdiev, Soloviev, Piersma and Klaassen2016). Conversely, changes in diet may result from a physiological response to abiotic climate that affects body size and foraging behavior (Pyke Reference Pyke1984; Dickman Reference Dickman1988). 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 Reference Levinton1983; Barnosky Reference Barnosky, Martin and Barnosky1993). 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. Reference Wood, Zelditch, Rountrey, Eiting, Sheets and Gingerich2007). 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. Reference Boyer, Puente, Gladman, Glynn, Mukherjee, Yapuncich and Daubechies2015). They allow for the possibility of discovering unanticipated patterns of shape change (Zelditch et al. Reference Zelditch, Swiderski and Sheets2012). 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. Reference Vitek, Manz, Gao, Bloch and Boyer2017). 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. Reference Adams, Rohlf and Slice2004). 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 Reference Kay and Hiiemae1974; Ungar et al. Reference Ungar, Healy, Karme, Teaford and Fortelius2018).
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 (P4–M3) 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 Reference Korth1979) 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. Reference Vitek, Manz, Gao, Bloch and Boyer2017). 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 Reference Bown and Schankler1982; Smith et al. Reference Smith, Bloch, Strait and Gingerich2002; Rose et al. Reference Rose, Chew, Dunn, Kraus, Fricke and Zack2012). 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 Reference Gingerich2001; Secord et al. Reference Secord, Gingerich, Smith, Clyde, Wilf and Singer2006). The beginning of the PETM is coincident with the start of the Wasatchian North American Land Mammal Age (Koch et al. Reference Koch, Zachos and Gingerich1992) and the Eocene Epoch (Aubry et al. Reference Aubry, Ouda, Duipuis, Berggen and Van Couvering2007). 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. Reference Bowen, Koch, Gingerich, Bains and Corfield2001; Gingerich Reference Gingerich2001; Gingerich and Smith Reference Gingerich and Smith2006). 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. Reference Rose, Chew, Dunn, Kraus, Fricke and Zack2012; Baczynski et al. Reference Baczynski, McInerney, Wing, Kraus, Bloch, Boyer, Secord, Morse and Fricke2013). 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. Reference Bowen, Koch, Gingerich, Bains and Corfield2001; Gingerich Reference Gingerich1989, Reference Gingerich2001). 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 Reference Gingerich2001). The Wa-1 biozone is preserved in only about 14 m of section in the southern BHB (Morse et al. Reference Morse, Chester, Boyer, Smith, Smith, Gigase and Bloch2019). To more fully characterize the post-PETM interval, we also include specimens from the subsequent Wa-2 biozone (Gingerich Reference Gingerich2001; 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. Reference Secord, Gingerich, Smith, Clyde, Wilf and Singer2006; 2, Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012; 3, van der Meulen et al. Reference van der Meulen, Gingerich, Lourens, Meijer, van Broekhuizen, van Ginneken and Abels2020; 4, Chew Reference Chew2009. Note that van der Meulen et al. Reference van der Meulen, Gingerich, Lourens, Meijer, van Broekhuizen, van Ginneken and Abels2020 measured a shorter total duration of the PETM than was used in Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012. 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. Reference Bowen, Bralower, Delaney, Dickens, Kelly, Koch, Kump, Meng, Sloan, Thomas, Wing and Zachos2006; Kraus et al. Reference Kraus, McInerney, Wing, Secord, Baczynski and Bloch2013). 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 M1s and M2s 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. Reference Boyer, Gunnell, Kaufman and McGeary2016; 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. Reference Rose, Chew, Dunn, Kraus, Fricke and Zack2012). 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. Reference Spradley, Pampush, Morse and Kay2017). 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.).
Sample sizes of the three focal taxa for each type of measurement. M1 sample sizes refer to specimens measured to calculate crown area. M2 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 M1 (Gingerich et al. Reference Gingerich, Smith and Rosenberg1982; Hlusko et al. Reference Hlusko, Lease and Mahaney2006; 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. Reference Yezerinac, Lougheed and Handford1992). 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. Reference Zelditch, Swiderski and Sheets2012), but it is less strongly correlated with body mass and therefore not preferred over log-transformed crown area in this study (Gingerich and Smith Reference Gingerich, Smith and Jungers1984; St. Clair and Boyer Reference St. Clair and Boyer2016). To test hypotheses of size change between different climatic or ecological regimes, we compared time bins using a bootstrapped comparison of M1 area, length, and width means with 1000 iterations (Kowalewski and Novack-Gottshall Reference Kowalewski and Novack-Gottshall2010). Based on previous morphometric studies, bins were not compared if either bin contained fewer than five specimens (Polly Reference Polly2003). Pairwise tests were corrected using the false discovery rate correction (Benjamini and Hochberg Reference Benjamini and Hochberg1995).
Molar Shape Measurement
The crown of M2 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. Reference Townsend, Richards and Hughes2003; Table 2). The complex occlusion of the interlocking molar cusps constrains M2 morphology to be relatively conservative due to the strong selective disadvantage of malocclusion in mammals with tall molar cusps (Gingerich and Winkler Reference Gingerich and Winkler1979).
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. Reference Evans, Wilson, Fortelius and Jernvall2007; Boyer Reference Boyer2008; Bunn et al. Reference Bunn, Boyer, Lipman, Clair, Jernvall and Daubechies2011; Pineda-Munoz et al. Reference Pineda-Munoz, Lazagabaster, Alroy and Evans2017; López-Torres et al. Reference López-Torres, Selig, Prufrock, Lin and Silcox2018). 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. Reference Evans, Wilson, Fortelius and Jernvall2007; Pineda-Munoz et al. Reference Pineda-Munoz, Lazagabaster, Alroy and Evans2017). 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 Reference Bown and Schankler1982; Smith et al. Reference Smith, Bloch, Strait and Gingerich2002; Beard and Dawson Reference Beard and Dawson2009; Rose et al. Reference Rose, Chew, Dunn, Kraus, Fricke and Zack2012). In contrast, RFI and DNE are sensitive to feature quality as well as feature quantity (King et al. Reference King, Arrigo-Nelson, Pochron, Semprebon, Godfrey, Wright and Jernvall2005; Gutzwiller and Hunter Reference Gutzwiller and Hunter2015; Pampush et al. Reference Pampush, Spradley, Morse, Griffith, Gladman, Gonzales and Kay2018). The RFI metric expresses crown height as a ratio of 3D occlusal area to 2D planimetric footprint area (Boyer Reference Boyer2008). The DNE metric integrates crown curvature across a digital surface, effectively capturing tooth sharpness (Bunn et al. Reference Bunn, Boyer, Lipman, Clair, Jernvall and Daubechies2011; Winchester et al. Reference Winchester, Boyer, Clair, Gosselin-Ildari, Cooke and Ledogar2014; Pampush et al. Reference Pampush, Spradley, Morse, Harrington, Allen, Boyer and Kay2016a). 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 Reference Evans and Sanson2003, Reference Evans and Sanson2005). 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 Reference Puente2013; Boyer et al. Reference Boyer, Puente, Gladman, Glynn, Mukherjee, Yapuncich and Daubechies2015). 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. Reference Pampush, Winchester, Morse, Vining, Boyer and Kay2016b) 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 M2 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. Reference Vitek, Manz, Gao, Bloch and Boyer2017). 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. Reference Yezerinac, Lougheed and Handford1992; Fruciano Reference Fruciano2016; Vitek et al. Reference Vitek, Manz, Gao, Bloch and Boyer2017).
M2 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. Reference Yezerinac, Lougheed and Handford1992).
Shape differences between time bins in the M2 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.)
Illustration of the univariate measurements taken post hoc on M2s. An M2 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. Reference Ungar, Healy, Karme, Teaford and Fortelius2018). 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 Reference Kay and Hiiemae1974). 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 Reference Green, Croft, Croft, Su and Simpson2018). 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. Reference Adams, Collyer, Kaliontzopoulou and Sherratt2017). A minimum sample size of N = 5 was imposed to test a bin for significant differences (Polly Reference Polly2003). 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 Reference Kowalewski and Novack-Gottshall2010). Pairwise tests were corrected using the false discovery rate correction (Benjamini and Hochberg Reference Benjamini and Hochberg1995).
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. Reference Smith, Bloch, Strait and Gingerich2002; Yans et al. Reference Yans, Strait, Smith, Dupuis, Steurbaut and Gingerich2006; Rose et al. Reference Rose, Chew, Dunn, Kraus, Fricke and Zack2012). 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 M1 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 M1 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 M1 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 M2 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 M2 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. Reference Dushoff, Kain and Bolker2019; 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).
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 (Reference Kay and Hiiemae1974), 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 M1 and crown surface topography of M2 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 Reference Clyde and Gingerich1998; Gingerich Reference Gingerich2006; McInerney and Wing Reference McInerney and Wing2011; Wing and Currano Reference Wing and Currano2013). 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. Reference Hunt, Bell and Travis2008). 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. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012). 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. Reference Bowen, Bralower, Delaney, Dickens, Kelly, Koch, Kump, Meng, Sloan, Thomas, Wing and Zachos2006). 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 Reference Kay and Hiiemae1974; Kay Reference Kay1975; Janis and Fortelius Reference Janis and Fortelius1988; Evans et al. Reference Evans, Wilson, Fortelius and Jernvall2007; Boyer Reference Boyer2008). 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 Reference Crompton1971; Kay and Hiiemae Reference Kay and Hiiemae1974). Greater amounts of wear or earlier appearance of wear facets implies that those regions of the tooth see greater use than others (Fortelius Reference Fortelius1985).
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 (Reference Kay and Hiiemae1974) and with shearing facets 2, 3, and 6 of the marsupial Didelphis marsupialis in the framework of Crompton and Hiiemae (Reference Crompton and Hiiemae1970). 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 Reference Kay and Hiiemae1974; Strait Reference Strait1993). 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 Reference Kay and Hiiemae1974). 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 Reference Strait1993). This correspondence supports the hypothesis that shape changes are localized on functionally important features and are therefore likely targets of selection (Barrett and Hoekstra Reference Barrett and Hoekstra2011). 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 Reference Loffredo and DeSantis2014; Green and Croft Reference Green, Croft, Croft, Su and Simpson2018).
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 M2 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 M1, size of the M2, morphology of the M3, 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. Reference Storz, Cheviron, McClelland and Scott2019). However, previous studies of molar morphology have found that it can closely track isotopically inferred changes in diet (Kimura et al. Reference Kimura, Jacobs and Flynn2013). 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 [Reference Clyde and Gingerich1994] and Wood et al. [Reference Wood, Zelditch, Rountrey, Eiting, Sheets and Gingerich2007]), 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 Reference Clyde and Gingerich1994; Wood et al. Reference Wood, Zelditch, Rountrey, Eiting, Sheets and Gingerich2007; Uyeda et al. Reference Uyeda, Hansen, Arnold and Pienaar2011). 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. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012). 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. Reference Smith, Betancourt and Brown1995) and equid body size during the PETM (Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012). 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 Reference Holt1990). 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 Reference Johansson2008). 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. Reference Churchfield, Nesterenko and Shvarts1999; Ochocinska and Taylor Reference Ochocinska and Taylor2003; Yom-Tov and Yom-Tov Reference Yom-Tov and Yom-Tov2005; White and Searle Reference White and Searle2007; but see Contoli et al. [Reference Contoli, Battisti and Buscemi2000] and Rychlik et al. [Reference Rychlik, Ramalhinho and Polly2006] for contrary expectations). Additional study of medium-sized herbivorous mammals, such as the “condylarths” Ectocion and Copecion (Gingerich Reference Gingerich2006), 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. Reference Smith, Betancourt and Brown1995; Secord et al. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012). 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. Reference Purvis, Gittleman, Cowlishaw and Mace2000; Cardillo and Bromham Reference Cardillo and Bromham2001; McCain and King Reference McCain and King2014). 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. Reference Secord, Bloch, Chester, Boyer, Wood, Wing, Kraus, McInerney and Krigbaum2012; Felibert et al. Reference Felibert, Morse, Strait, Boyer and Bloch2017). 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 Reference Meiri and Dayan2003; Smith et al. Reference Smith, Brown, Haskell, Lyons, Alroy, Charnov and Dayan2004). 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 Reference Schmidt-Nielsen1984; Terry and Rowe Reference Terry and Rowe2015). 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. Reference Purvis, Gittleman, Cowlishaw and Mace2000).
Acknowledgments
Thanks to J. Bourque and A. Poyer for specimen prep; T. Bown for locality and meter-level data; M. Clementz, A. Millhouse, N. Pyenson, K. Rose, and L. Vietti, for specimen access; P. Holroyd and R. Hulbert for cataloging assistance; A. Wood and T. Gao for coding assistance; G. Schieffle, J. Thostenson, and E. Stanley for CT-scanning assistance at Duke University and the University of Florida; K. Rosenbach for assistance segmenting CT scans; six anonymous reviewers for feedback that improved the article; and many iterations of the Bighorn Basin field crew for specimen-collection assistance. This research was supported by National Science Foundation grants DGE-1315138 (N.S.V.), BCS-1440558 (J.I.B.) BCS-1440742 (D.M.B.), EAR-0640076 (J.I.B., J. Krigbaum, R. Secord), EAR-0719941 (J.I.B.), DEB-0208281 (S.G.S.), and BCS-1552848 (D.M.B.). Vertebrate fossils were collected under the Bureau of Land Management permits to J.I.B. (PA04-WY-113, PA10-WY-185). This is University of Florida Contribution to Paleobiology 857.
Data Availability Statement
Data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.mkkwh70zr; and the GitHub Digital Repository: https://github.com/nsvitek/PETM-erinaceomorphs.
Supplementary Figure 1. Cropping repeatability for each species: Colpocherus (A, B), Macrocranion (C, D), and Talpavoides (E, F). Location of replicate surfaces on PC 1 and PC 2 in comparison to overall dataset (A, C, E). Distribution of repeatability across PCs (B, D, F).
Supplementary Figure 2. Patterns of M1 length (A−C) and width (D−F) through time for Colpocherus (A, D), Macrocranion (B, E), and Talpavoides (C, F).
Supplementary Figure 3. The first two PCs resulting from a PCA of 2048 pseudolandmarks spread evenly across the M2 samples for each species: A, Colpocherus, B, Macrocranion, and C, Talpavoides. Plotting character shapes and colors indicate sampling bins.
Supplementary Figure 4. Dietary proxies for the M2 among sampling bins for each species. Relief index (RFI; A−C) and Dirichlet normal energy (DNE; D−F) in Colpocherus (A, D), Macrocranion (B, E), and Talpavoides (C, F). Dashed lines indicate the beginning and end of the carbon isotope excursion delimiting the PETM. Plotting character shapes and colors indicate smallest sampling bins. Light brackets indicate bins being compared. Heavier brackets indicate statistical comparisons (not performed for N < 5 within a bin; none significantly different).
Supplementary Fig. 5. Examples of increasing tooth wear in Macrocranion (A−C) and Talpavoides (D−F). A, UF 283308, B, UF 283317, C, UF 283334, D, UF 330408, E, UF 434082, and F, UW 9688. Scale bars at the bottom of each column equal 1 mm. 1, wear facet one; 3, wear facet three; 4, wear facet four of Kay and Hiiemae (Reference Kay and Hiiemae1974); a, anterior; b, buccal; co, cristid obliqua; ent, entoconid; hyp, hypoconid; hypc, hypoconulid; l, lingual; met, metaconid; p, posterior; paco, paraconid; postcr, postcristid; proco, protoconid; procr, protocristid; talba, talonid basin; taln, talonid notch; triba, trigonid basin.
Open access
