INTRODUCTION
Small mammal communities from the last glacial period are composed of the same species found in modern environments; however, communities and individual ranges of species have undergone dramatic changes in response to past environmental change (Graham et al., Reference Graham, Lundelius, Graham, Schroeder, Toomey, Anderson and Barnosky1996; Grayson, Reference Grayson2006; Blois et al., Reference Blois, McGuire and Hadly2010). Many possess unique, nonanalog associations of taxa whose modern geographic ranges do not overlap (Graham and Lundelius, Reference Graham and Lundelius1984). Late Quaternary small mammal faunas document ecological change and biotic responses to past and novel climates, and their study informs our understanding of community responses to environmental change more broadly (Blois et al., Reference Blois, Zarnetske, Fitzpatrick and Finnegan2013). Although North America is one of the most heavily studied regions for late Quaternary paleoecology, well-dated, temporally long vertebrate faunal sequences from the northern Great Plains of North America are rare (Graham, Reference Graham1987; FAUNMAP Working Group, 1994).
Here, we present analyses of small mammal (<2 kg) community composition from Don’s Gooseberry Pit (DGP), a cave in the southeastern part of the Black Hills of South Dakota (Fig. 1). DGP contains an approximately 18,000-yr-old faunal sequence in this data-deficient region. We assess issues of time averaging within the stratigraphic sequence and present analyses of small mammal community turnover. We ask if nonanalog faunas are present and provide high-quality radiocarbon dates to constrain the timing of occurrence for some ecologically and climatically significant taxa.
Location of Don’s Gooseberry Pit (DGP) in the northern Great Plains and other fossil sites of Dicrostonyx in North America. The site is found in the southeast corner of the Black Hills on Forest Service land. The modern distribution of lemming species of Dicrostonyx (shaded, stippled, and cross-hatched areas) are shown (Patterson et al., Reference Patterson, Ceballos, Sechrest, Tognelli, Brooks, Luna, Ortega, Salazar and Young2003). Ancient DNA (aDNA) for specimens from DGP and Prairieburg suggest that all fossils south of the ice sheet are probably related to the modern species D. richardsoni (gray shaded range).
Regional setting and description of the site
The Black Hills in southwestern South Dakota and eastern Wyoming were situated south of the terminal limits of the Laurentide Ice Sheet during the last glacial period. Their higher elevation (~2182 m) relative to the surrounding plains (~910 m) creates an isolated “island” of forest surrounded by a “sea” of grasslands. The crystalline rock core of the Black Hills is surrounded by concentric exposures of Paleozoic, Mesozoic, and Cenozoic rock. The Mississippian-age Paha Sapa Limestone is extensively karstified. Many of these caves contain sediments that preserve vertebrate fossils, including DGP, which is located in the southeastern Black Hills about 20.5 km north of Hot Springs, South Dakota. The cave is on a relatively flat ridge top that trends north. Steep-sided valleys that originate near the cave are on the east and west sides of this ridge. The cave is surrounded by Ponderosa pine (Pinus ponderosa) forest interspersed with open grassy areas and lies at roughly 1667 m elevation.
The modern vegetation types of the Black Hills form an elevational gradient. Boreal taxa like spruce (Picea glauca) at the higher elevations intermix with Ponderosa pine forests at midaltitudes (Hoffman and Alexander, Reference Hoffman and Alexander1987). Open grassy areas, some quite large like Reynold’s Prairie, are interspersed with the forest. At lower elevations, deciduous taxa like oaks (Quercus spp.) intermingle with, and eventually give way to, grasslands. Various mammal species are associated with, and indicative of, these different vegetation types (Turner, Reference Turner1974).
The distribution of vegetation fluctuated in and around the Black Hills with the climate oscillations of the late Quaternary. Broader expanses, and greater connectivity, of spruce woodland across the northern Great Plains during the Pleistocene (Watts and Wright, Reference Watts and Wright1966; Rhodes, Reference Rhodes1984; Stewart, Reference Stewart1987) facilitated the interchange of forest-adapted plants and animals from the Black Hills, Rocky Mountains, and boreal environments farther north. As the environment became warmer and drier throughout the Holocene, grasslands expanded and forested habitats retreated to higher elevations. Forest-adapted species then became isolated in a Black Hills refugium “island,” as they are today. The fauna from DGP documents these climate and environmental fluctuations.
MATERIALS AND METHODS
Treatment of specimens
Sediment from DGP was excavated by trowels in 10-cm and 5-cm levels within natural stratigraphic units and then wet sieved to recover fossils (see the Supplementary Materials for details of excavation and processing techniques). Identifications were based solely on morphological characters, primarily of dentition (see the Supplementary Materials for details), and not the proximity of DGP to modern geographic ranges. This avoids the circular logic of using modern biogeography to inform the identity of a specimen that will then be used to constrain past distributions (Bell et al., Reference Bell, Gauthier and Bever2010). When possible, genera were placed into morphological “clades” of candidate species (Supplementary Materials). In the case of Dicrostonyx, ancient DNA (aDNA) identified individuals as Dicrostonyx richardsoni (Richardson’s lemming) from DGP as well as other late Pleistocene sites (Fig. 1; Fulton et al., Reference Fulton, Norris, Graham, Semken and Shapiro2013). The number of identified specimens (NISP) for each taxon per sample was tabulated (Supplementary Table 1); these data have been made available through the Neotoma Paleoecology Database (http://www.neotomadb.org). Samples were designated by excavation level and stratigraphic unit. Minimum number of individuals (MNI) is conventionally used as the counting unit in most paleoecological analyses. We are deviating from this convention here, as MNI is highly dependent on the number of subdivisions in an excavation site (in this case, by level and sedimentary units). As excavated subdivisions become finer, MNI approaches the value of NISP total for the entire site (Grayson, Reference Grayson2014).
Cluster analysis and community turnover
Our statistical analyses were conducted on rodents, lagomorphs, and shrews (Sorex sp.), although the complete mammal fauna also includes artiodactyls and carnivores in low numbers (Table 1). Counts were square-root transformed to emphasize the presence of rare taxa, and then the small mammal fauna was clustered in two ways using agglomerative hierarchical clustering (Maechler et al., Reference Maechler, Rousseeuw, Struyf, Hubert and Hornik2011) in the programming language R (R Development Core Team, 2009). Samples were first clustered by species composition using a Euclidean distance metric to find natural groupings of taxa to produce “habitat groups.” The data were then transposed and clustered, using the same algorithm and distance metric, to identify “species functional groups” clustered together by samples.
Complete mammal fauna from Don’s Gooseberry Pit as of 2014. Statistical analyses were conducted on the rodents, shrews, and lagomorphs. Details on number of identified specimens and their stratigraphic location are provided in the Supplementary Materials.
We assessed community turnover in two ways: first, within the context of our interpretation of the cave stratigraphy and then by directly radiocarbon dating ecologically diagnostic taxa from throughout the excavation. For the first assessment, the percent composition of taxa in each species functional group was plotted over the stratigraphic sequence. Percent composition was calculated by dividing the NISP found in species functional group by the total NISP for the sample. For the second assessment, we radiocarbon dated three taxa (D. richardsoni, Microtus with a lower first molar (m1) morphology possessing three closed reentrant angles, and Myodes sp.) using purified bone collagen and calibrated the dates to IntCal13 (Reimer, Reference Reimer2013). The dates were obtained from specimens with diagnostic teeth, so an exact date could be associated to a specific taxon.
RESULTS
Stratigraphy, chronology, and time averaging
Supplementary Table 2 describes the sedimentary units from DGP. Most of the sediment was washed into DGP on its northeastern edge where a waterworn groove is present on the limestone lip of the cave entrance. In general, the sediments in the various recognized stratigraphic units are a heterogeneous and mostly an unsorted mixture of muds, sands, gravels, angular cobbles, and boulders. Although we have not conducted particle-size analyses of the deposits, we did recognize and record, during excavation, varying sizes of the larger clasts (cobbles and boulders) in various stratigraphic units.
None of the stratigraphic units appear to dip significantly as is common in debris cones (Supplementary Figs. 1 and 2). We believe the small size of the cave constrained development of a large angular cone. The horizontal nature of the bedding is confirmed by pebble zones that can be traced along bedding planes without showing any breaks or erosion features (Supplementary Figs. 2, 3, and 4). There are some lens-shaped cut-and-fill deposits (Supplementary Fig. 2), but each natural stratigraphic unit was kept separate (Supplementary Materials). There are a few rodent burrows, but they are distinct in sediment color and their consistency (Supplementary Fig. 3). Burrows were excavated separately, removed, and not included in analyses.
There is considerable time averaging in the fauna of DGP. Radiocarbon dates (methods of dating are discussed in the Supplementary Materials) from specific taxa range from the full glacial period (18,235 cal yr BP) through the middle Holocene (6287 cal yr BP) (Table 2). Examination of the dates and the taxa demonstrate that they are temporally mixed within stratigraphic units and excavation levels (Fig. 2). The degree of time averaging appears to be extensive in level 10, unit 4, as more than 5000 yr are represented in this sample. The temporal mixing precludes calculating an explicit age-depth model for DGP, although there is a clear sequence of replacement (Dicrostonyx-Myodes-Microtus) when the dates are arranged chronologically in sequence (Fig. 3). When superimposed with δ18O values from the Greenland ice core (Johnsen et al., Reference Johnsen, Clausen, Dansgaard, Gundestrup, Hammer, Andersen and Andersen1997), this sequence follows trends in warming temperature (Fig. 3), with cold-adapted taxa such as D. richardsoni restricted to cold climates of the full and late glacial periods, and warmer taxa such as Microtus restricted to the Holocene.
Plot of radiocarbon dates from Don’s Gooseberry Pit by taxon, level, and age. Although three taxa occur within levels between 133 and 143 cm, the dates for all three are diachronous, indicating that the fauna is temporally mixed. The same pattern is apparent for the levels below and above this one. Age error is obscured by data point markers but is provided in Table 2.
Plot of median calibrated dates from individual species. Calibrated ages are given along with δ18O values from the Greenland Ice Core Project (GRIP) ice core in Greenland as a proxy for climate (Johnsen et al., Reference Johnsen, Clausen, Dansgaard, Gundestrup, Hammer, Andersen and Andersen1997; https://www.ncdc.noaa.gov/paleo/study/17839 [accessed September 20, 2017]). Colder temperatures are indicated by lower δ18O values. Age error is obscured by data point markers but is provided in Table 2. Oxygen isotope values were interpolated from the ice core data for the dates of each taxon to superimpose their position on the δ18O record. The approximate duration of the Younger Dryas from 11,500 to 12,700 cal yr BP is shown in black.
Radiocarbon dates from Don’s Gooseberry Pit. Dates were run on purified collagen from jaws containing taxonomically diagnostic teeth. Dates were calibrated to IntCal13 (Reimer, Reference Reimer2013).
Community composition and ecological turnover
Clustering of samples revealed structure throughout the excavation and produced four distinct groupings (Fig. 4). Likewise, clustering of taxa produced five clearly defined functional groups (Fig. 5). Distinct faunal communities appear throughout the excavation. Cluster 1 contains “Ubiquitous” taxa that occur throughout most samples (Supplementary Table 1). Cluster 2 contains only red-backed voles (“Myodes”). Cluster 3 contains two taxa: Microtus belonging to a morphological clade that has five to seven triangles in the lower first molar and Dicrostonyx. Of these, Dicrostonyx is more abundant; we are therefore referring to this grouping as the “Tundra” cluster. Cluster 4 contains the greatest diversity of taxa and can be further subdivided into a fauna dominated by squirrels (Sciurus sp.) and other taxa that may be found in closed environments (4a, “Forest”) and a fauna that contains a mixture of species found in a variety of environments today (4b, “Mixed”; e.g., forested and open environments). Cluster 4b contains many taxa that are relatively rare throughout most of the strata (Supplementary Table 1). Cluster 5 is characterized by taxa that may be found in dry or arid environments, although there is some potential variation in habitat preference within genera; this cluster is called the “Arid” species cluster.
Sample cluster analysis. Clusters are labeled and alternately highlighted in gray. Natural groupings occur because of similarity of shared taxa.
Species cluster analysis. Natural groupings occur from species that are associated together in samples. Species functional groups are labeled and alternately highlighted in gray. Group 1 consists of “Ubiquitous” taxa, group 2 consists entirely of Myodes sp., group 3 consists of glacial taxa, group 4a is the “Forest” group, group 4b is the “Mixed” group, and group 5 is the “Arid” group.
Aligning the sample clusters (Fig. 4) and the species clusters (Fig. 5) with percent composition reveals the taxonomic associations that are driving community structure throughout the excavation (Fig. 6). The greatest diversity of species in habitat cluster A1 is from the 4a, “Forest,” species cluster, as well as a large percentage of Onychomys from the “Arid” cluster 5. Cluster A2 is primarily characterized by “Ubiquitous” taxa, Myodes, and Dicrostonyx. Cluster B1 consists of Microtus, Myodes, Dicrostonyx, and Synaptomys cooperi, whereas the diversity of species occupying cluster B2 constitutes a combination of “Ubiquitous” taxa and Myodes, plus the “Mixed” cluster (4b) and “Arid” taxa (cluster 5). Clusters B3i and B3ii are similar; however, B3ii lacks Myodes, and B3i possesses a greater diversity of species from the “Mixed” species cluster. Cluster C is composed of a great diversity of species from across all species functional groupings. Cluster D is also quite diverse but lacks the density of Dicrostonyx that is present in cluster C.
The relationship between habitat groups and species functional groups. Dot sizes are proportional to the relative percentage of each taxon in a sample.
Stratigraphic ordering of species functional groups across samples reveals an overall warming and drying environmental trend (Fig. 7). The fauna from the lowest levels (L14-U8) are thought to represent full glacial age deposits and are composed of approximately 48% of “Tundra” taxa. A shift occurs at level 13 where Myodes becomes dominant at more than 70% of diversity. The community changes again at level 12, unit 2, where “Forest” and “Arid” taxa dominate, at 60% and 35%, respectively. Between levels 12 and 10, community composition switched to one composed mostly of Myodes (24%) and “Ubiquitous” taxa (59%). Younger levels are dominated by “Ubiquitous” (>50%) and “Arid” (>10%) taxa that are found around the cave today.
Changes in the percentage composition of the Don’s Gooseberry Pit (DGP) fauna by species functional groups throughout excavation units. Samples contain 10 or more species. Samples obtained above level 10 are composed primarily of “Ubiquitous” and “Arid” taxa (see Fig. 4). Samples obtained below level 10 are dominated by Myodes and “Ubiquitous” and “Tundra” species. Level 10 in the DGP excavation likely represents the Holocene-Pleistocene boundary.
DISCUSSION
Stratigraphy, chronology, and time averaging
The source of the time averaging in DGP is unclear; however, the apparent horizontal nature of the deposits in DGP suggested that they had not been overturned by extensive burrowing. Cut-and-fill deposits always have the potential for mixing, and although these types of deposits occurred, they were not dominant. Excavation techniques have been careful to separate natural stratigraphic units and intrusive structures like rodent burrows to minimize the potential for temporal mixing (Supplementary Materials).
Initially, we had inferred from their appearance that the deposits and fauna were contemporaneous. Mixing was only apparent after differing radiocarbon dates were obtained from multiple taxa within the same stratigraphic levels and units (Fig. 2). DGP serves as a cautionary note about accepting the age of a fauna on the basis of a single date or a series of dates on a specific taxon. Single dates from a stratigraphic horizon will not capture the potential time averaging of the stratum. Likewise, a series of dates on a specific taxon that may have occupied the area for a short interval of time (e.g., Dicrostonyx) will not capture the diachronic mixture of the entire fauna.
Community composition and ecological turnover
Many species in the DGP fauna (Table 1) are extralimital, which is characteristic of most late Pleistocene faunas (Graham and Lundelius, Reference Graham and Lundelius1984). Their presence reinforces past tundralike or tundra-boreal forest ecotone paleoenvironmental reconstructions for the Black Hills. Examples include pika (Ochotona sp.), which is primarily found in alpine tundra or at higher elevations in mountainous regions, and the boreal bog lemming (Synaptomys borealis) that lives in meadows in spruce forests farther north or in the Rocky Mountains (Banfield, Reference Banfield1974; Turner, Reference Turner1974; Higgins et al., Reference Higgins, Stukel, Goulet and Backlund2002). Myodes is also cold adapted and can be found in the Black Hills today, but only at higher elevations (Turner, Reference Turner1974; Graham, R.W., unpublished trapping records 2006–2009). Microtus from the five to seven closed triangle m1 morphoclade (Supplementary Materials) could represent several taxa found in a variety of environments including tundra, spruce forest, and grassland.
The recovery of Dicrostonyx from DGP deserves special notice. DGP is the first record of this genus for the state of South Dakota (Fig. 1), and ours are the first direct radiocarbon dates for Dicrostonyx (Fulton et al., Reference Fulton, Norris, Graham, Semken and Shapiro2013); previous studies have relied on associated dates (Guilday et al., Reference Guilday, Martin and McCrady1964; Zeimens and Walker, Reference Zeimens and Walker1974; Grady and Garton, Reference Grady and Garton1981; Foley, Reference Foley1984; Voorhies and Corner, Reference Voorhies and Corner1985; Bonnichsen and Oliver, Reference Bonnichsen and Oliver1986; Chomko and Gilbert, Reference Chomko and Gilbert1987; Foley and Raue, Reference Foley and Raue1987; Graham, Reference Graham1987; Miller, Reference Miller1987; Semken and Falk, Reference Semken and Falk1987; Walker, Reference Walker1987). aDNA of fossil Dicrostonyx obtained from DGP and Prairieburg (Iowa, USA) indicate those specimens are most closely related to a modern species, D. richardsoni (Richardson’s lemming) (Fulton et al., Reference Fulton, Norris, Graham, Semken and Shapiro2013) (Fig. 1). Dicrostonyx is found distributed along the front of the Laurentide Ice Sheet in other states (Wyoming, Nebraska, Iowa, Wisconsin; Fig. 1) during the full glacial period (Mead and Mead, Reference Mead and Mead1989), though most of these sites have fewer than 10 specimens. At DGP, we have found dozens of specimens within multiple stratigraphic units, which raises an interesting, but as yet unanswered, question as to why they are so abundant at this cave.
Biogeographically, the wide late Pleistocene geographic distribution of D. richardsoni stands in stark contrast to its more relict distribution along the western edge of Hudson Bay today and the wide arctic distribution of other Dicrostonyx species (Fig. 1). Furthermore, this reduction in geographic range suggests that D. richardsoni may be vulnerable to extinction with future global warming. A refugium of these animals could have remained in the Black Hills longer than at lower elevations and latitudes. Because we did not conduct aDNA on all specimens, other specimens of Dicrostonyx from DGP might be other species. However, because aDNA studies place the distribution of D. richardsoni at several different sites widely along the front of the Laurentide Ice Sheet during the late Pleistocene (Fig. 1; Fulton et al., Reference Fulton, Norris, Graham, Semken and Shapiro2013), we believe this is unlikely.
Despite the significant time averaging, the faunal composition is structured (Figs. 4–7) and changes in DGP are consistent with other late Quaternary mammal faunas (e.g., Graham, Reference Graham1987) and vegetation records (Valero-Garcés et al., Reference Valero-Garcés, Laird, Fritz, Kelts, Ito and Grimm1997; Grimm, Reference Grimm2001; Clark et al., Reference Clark, Grimm, Donovan, Fritz, Engstrom and Almendinger2002) from the northern plains. The fauna from the lowest excavation levels represents a colder environment and is dominated by species from the “Tundra” and “Myodes” clusters (Fig. 7). “Tundra” gives way to a fauna dominated by “Myodes,” which indicates it was somewhat warmer, but still a closed environment (Fig. 7). A major change in the fauna occurs in level 12; “Tundra” and “Myodes” are replaced by taxa from the “Forest” and “Arid” cluster, which may represent a mosaic of forest and grassland mixed habitats (Fig. 7). The fauna changes dramatically after level 10, becoming dominated by “Ubiquitous,” “Mixed,” and “Arid” species, suggesting an overall warming and drying trend (Fig. 7). These changes are consistent with an overall trend of late Quaternary warming (Fig. 3).
Assessment of nonanalog fauna and ecological replacement
Strategic radiocarbon dating indicates the dated taxa from DGP are diachronous (Figs. 2 and 3), and we find no evidence to support nonanalog associations. Although radiocarbon dating of individual taxa at other sites has documented biologically meaningful nonanalogs (Stafford et al., Reference Stafford, Semken, Graham, Klippel, Markova, Smirnov and Southon1999; Semken et al., Reference Semken, Graham and Stafford2010), DGP clearly demonstrates that for a fauna to be deemed contemporaneous, the taxa must be directly dated. It is possible that future dating of specimens from DGP may reveal contemporaneous nonanalogs.
Temporal ordering by radiocarbon dates of individual species documents the replacement of ecologically significant small mammal species that is consistent with the Greenland ice core record (Fig. 3; Johnsen et al., Reference Johnsen, Clausen, Dansgaard, Gundestrup, Hammer, Andersen and Andersen1997). Dicrostonyx specimens dating 17,000–14,000 cal yr BP are all from the full and late glacial periods. Specimens of Myodes date to the late glacial period (ca. 13,000 cal yr BP) but are younger than any of the Dicrostonyx specimens. Myodes does not occur in the Holocene at DGP. Throughout the early and middle Holocene (ca. 8 ka), as reflected in the appearance and persistence of the three closed triangle clade of Microtus, the climate became warmer. The current pattern of dates and taxa (Fig. 3) documents ecological replacement. This process can be, and probably is, the result of an individualistic species response, as species move in and out of the community (Figs. 6 and 7).
CONCLUSIONS
An overall pattern of environmental change at DGP is preserved, even though any particular stratum or excavation unit may not necessarily contain strictly contemporaneous taxa. For broadscale reconstructions beyond the age of radiocarbon dating, such time averaging may be tolerated, but researchers must not overextrapolate from their findings. DGP serves as a warning that it may not be possible to verify nonanalog biotas beyond radiocarbon dating limits (e.g., 40 14C ka BP). For fine-scale reconstructions like those of the late Quaternary, time averaging of more than a few hundred years may not be acceptable without extensive radiocarbon dating.
Individually dated specimens from DGP document temporal turnover of paleoecologically sensitive species and ecological replacement (Dicrostonyx-Myodes-Microtus) through time. As revealed by individually dated taxa, this replacement reflects the initial cold stage of the full glacial period, which is then replaced with continual warming throughout the late glacial period and early and middle Holocene. Radiocarbon dates on individual specimens of Dicrostonyx from DGP suggest that this taxon appears to have survived late in the Black Hills (17,083 cal yr BP), which may have served as an altitudinal refugium for D. richardsoni and other cold-adapted taxa.
Acknowledgments
The U.S. Forest Service (USFS) is thanked for its support and allowing the research project to occur on USFS property. Megan Buehler is especially thanked for helping us obtain our permits for excavation and travel across forest land. Personnel from the USFS were also instrumental in training us for sanitizing our field gear and equipment against white nose syndrome contamination. These efforts are appreciated. We wish to thank all of the volunteers from the Western Interior Paleontological Society for their efforts in providing labor, equipment, funds, and fun for the project. Without their help, this project could not have been done. We are grateful to Don Brandenburg and Angela Mathias who conducted all of the excavations in DGP. Special thanks to Eric Grimm who has served as a co–principal investigator on the project. We also thank volunteers from the Illinois State Museum and Penn State University (PSU). The Earth and Mineral Sciences Museum at PSU and the Illinois State Museum provided time for Russell Graham and Eric Grimm, respectively, to participate in the project. They also supplied some equipment and supplies. We are especially indebted to all of the PSU students who have worked on the project. We thank Doug Backlund from the South Dakota Department of Game, Fish and Parks for permission to trap small vertebrates in the vicinity of the cave to better sample the modern biota. We also want to thank the Guilday Hall volunteers for their time and effort in picking, literally, tons of screen wash matrix. Their dedication to the project is simply commendable. We thank three anonymous reviewers for their detailed comments, Curtis Marean and Derek Booth for helpful editorial comments, and Pat Shipman for proofreading the final manuscript. However, we assume all responsibility for interpretations and any errors. National Science Foundation grant EAR 0948652 supported MIP.
Supplementary material
For supplementary materials referred to in this article, please visit https://doi.org/10.1017/qua.2017.100
