Isotopic paleoecology of Northern Great Plains bison during the Holocene

Bison (Bison bison) are one of the few terrestrial megafauna to survive the transition into the Holocene and provide a unique opportunity to study a species on a broad spatiotemporal scale. Today, bison are primarily managed in small and isolated herds with little known about their ancestral ecology. We studied the carbon and nitrogen isotopes of Northern Great Plains bison from the terminal Pleistocene and throughout the Holocene to gain insight into their paleoecology. This time span is contemporary with the first population bottleneck experienced by bison at the end of the Pleistocene and includes the second bottleneck which occurred in the late 19th century. Results were compared with modern bison herd isotopic values from Theodore Roosevelt National Park (TRNP). Patterns of isotopic variation found in bison over time indicate significant (δ13C p = 0.0008, δ15N p = 0.002) differences in diet composition and correlate with climate throughout the Holocene. Isotopic relationships described here reveal the plasticity of ancient bison in unrestricted rangelands during periods of climatic fluctuations. Managers at TRNP and elsewhere should pursue opportunities to expand bison range to maximize forage opportunities for the species in the face of future environmental change.

help shed light on bison diet and foraging habitat selection over time, providing insights regarding physiological plasticity of the species relevant to management.
Isotopic biogeochemistry of collagen found in ancient bones and teeth is increasingly used in the construction of paleoecology and paleoenvironments. In addition to recording climatic variables such as temperature and precipitation, stable isotopic signatures encapsulate feeding strategies of animals from the past [11][12][13] . Isotopic ratios of carbon (δ 13 C) and nitrogen (δ 15 N) are assimilated into herbivore skeletal collagen and tooth dentin through diet, recording the isotopic composition of plant material consumed 14 . These values will remain the same over time in well-preserved specimens 15 . Bone collagen has a slow turnover rate and will record several years of the animal's life [16][17][18] . Tooth dentin forms at a specific point in development and records a shorter window of time, in the order of months 19,20 . The carbon isotopic signature depends on the proportion of plants having C 3 and C 4 photosynthetic pathways in ecosystems 21,22 as well as changes in atmospheric CO 2 from canopy cover in forested areas 23 or post-industrial revolution CO 2 emissions 13 . The values of C 3 plants in North America range from −30 to −24‰ and C 4 plants fall between −15 and −11‰ 24 . This allows us to distinguish between grazers, browsers, and mixed feeders due to the depletion in 13 C (i.e. lower δ 13 C values) observed in grazing diets 17 . Bison exhibit an enrichment factor of 6.3‰ 25 when carbon isotopes are assimilated into their skeletal tissue and this value needs to be factored into calculating the percentage of C 3 and C 4 plants in bison diet. Therefore, the expected δ 13 C values of bison feeding primarily on C 3 plants would be between −23.7 and −17.7‰, while C 4 bison diets would range from −8.7 and −4.7‰. Throughout the Holocene, the Great Plains have been largely dominated by Poaceae (grass) communities 26 , which can exhibit C 3 or C 4 photosynthetic pathways. The climate in the Northern Great Plains has predominately favored the C 3 subfamily Pooideae with a smaller amount of C 4 subfamilies, Panicoideae and Chloridoideae. The abundance of C 4 grasses increases as warm seasons get longer, allowing us to capture climatic changes in the diets of grazers 22 . Nitrogen isotopic values provide insight into moisture level and nutritional stress due to an observed increase in δ 15 N in animal tissue from the recycling of urea under conditions of drought 15,17,27,28 . However, there are other factors that contribute to nitrogen values in herbivores. Higher nitrogen can indicate warmer temperatures and a diet composed of more graminoids and herbs than trees and shrubs 13,[29][30][31][32] . In modern European bison, it was found that canopy cover had the biggest influence on δ 15 N, where less light will decrease nitrogen values in plants 29 .
North American isotopic studies of bison to date have primarily focused on Pleistocene paleoecology 30,31,[33][34][35][36] , climatic interpretations 11,32 , and values from modern herds 12,27,37 . At present, there are few isotopic studies of Holocene bison in North America. Existing research covers relatively short time periods or small sample sizes 11,[38][39][40][41] . The limited data on North American bison from the Holocene may be in part due to the prior perception of a relatively stable climate during this epoch, though the most recent studies of paleoclimate portray the Holocene as a dynamic period with fluctuations in temperature and precipitation [42][43][44] . An analysis of approximately fifty paleoclimatic records of greenhouse gases, glacial coverage, and pollen profiles determined that the Holocene had several periods of sudden climate change outlined by variations in atmospheric circulation, moisture, and temperature changes 43 . For instance, a study of pollen records found three periods in which climate rapidly heated during the Holocene 43 and evidence of a steep drop in temperature at ~8.2 thousand years ago as well as other large fluctuations in temperature throughout the Holocene are supported by the analysis of Greenland ice cores 44 . Therefore, bison in North America during the Holocene had to adapt to a wider range of climatic conditions than previously thought. The amount of variation seen in bison isotopic values throughout specific time periods provides insight into their use of different resources as the environment changed 39,40 .
The above-described wealth of information available on Holocene and Pleistocene isotope variation provides a unique opportunity to better understand historic bison through analysis of their bones and teeth represented in archaeological and natural history collections in the Northern Great Plains. Here, we explore a subset of that history through analysis of bison specimens from 22 archaeological sites across the Northern Great Plains spanning the Late Pleistocene to the Late Holocene in comparison to extant bison from Theodore Roosevelt National Park (TRNP), North Dakota, USA (Fig. 1). Our objectives are to: 1) identify δ 13 C and δ 15 N variation indicative of environmental change across time, and 2) utilize isotopic signatures to elucidate feeding ecology of historic bison in context of modern counterparts in North Dakota. We expect that isotopic values in bison will follow changes in climate and subsequently, isotope values of plant material. We expect increased δ 13 C and δ 15 N during periods of warm and dry conditions 22 .

Results
Fifty-five of the fifty-nine samples of ancient bison bones and teeth yielded sufficient collagen (>1%) to produce reliable data. All samples included in the analysis have C:N ratios between 2.9 and 3.6, relative carbon (%C) >30%, and relative nitrogen (%N) between 11 and 16%, indicative of well-preserved collagen (Table 1) 45,46 . Five modern bison tooth samples (Bison 100-Bison 104) from TRNP herds returned carbon and nitrogen isotopic ratios for comparison with ancient samples.
As previously stated, dentin and bone collagen represent different lengths of time in the bison's life. The sampling method in this study for dentin likely captures less than a year 20 while the collagen from bone samples describes the average over several years of the animal's diet [16][17][18] . To compare the variation found in each tissue type, we conducted t-tests for the dentin and bone samples as a whole and also separated into temporal groups. We found no significant differences in the means or variation found in sample tissue types (Supplementary Figure S1), therefore, they are treated functionally the same for the purpose of this paper.
Modern bison had the least variation in carbon and nitrogen values and the lowest mean δ 15 N value, 4.9‰. The highest variation for both isotopes was seen in Middle Holocene bone samples. The highest mean δ 13 C and δ 15 N values, −16.8‰ and 8.2‰ respectively, were observed in Middle Holocene bison dentin. Late Pleistocene and modern bison dentin were the most depleted in 13 C with a mean δ 13 C of −20.6‰ and −20.5‰, respectively (  (Table 2). Late Pleistocene bison consumed entirely C 3 vegetation with an increase in the consumption of C 4 plants in the Early Holocene, a peak abundance of C 4 material in the Middle Holocene and a small decrease again in the Late Holocene (Table 2, Fig. 2a). Modern bison in TRNP exhibited an entirely C 3 diet, congruent with the vegetation found in the park 48 .

Discussion
This study uses the archaeological record of bison in the Northern Plains to understand their evolutionary responses to environmental change and provide insight for best practices in bison conservation and management. Our analysis of isotopic values from bison remains spans the terminal phase of the Late Pleistocene to present, representing post-glacial changes in bison diet and vegetation associated with changing climate during the recent natural history of the species. The sample assemblage is contemporary with two bison population bottlenecks. The first bottleneck occurred during the terminal Pleistocene and the second in the late 19 th century when bison were nearly extirpated by humans 6,9 . The Late Pleistocene in the Northern Great Plains is described as a time of sudden environmental change and a significantly wetter landscape after the recent retreat of the Laurentide ice sheet 49 . Pollen records from the study region indicate an abundance of Picea (evergreen) species and relatively low amounts of herbaceous plants (grasses and forbs) 49,50 . Late Pleistocene bison samples in the current study are derived from Beacon Island, a Paleoindian kill site in the Agate basin (Fig. 1). The concurrent stratigraphic layer at this site exhibits C 3 dominated plant material 51 and is congruent with the 100% C 3 diet recorded in the δ 13 C of bison remains ( Table 2).  www.nature.com/scientificreports www.nature.com/scientificreports/ However, bison δ 15 N values from the Late Pleistocene appear surprisingly low considering environmental conditions in the Late Pleistocene were adverse enough to wipe out the majority of megafaunal species and facilitate bison's first recorded population bottleneck 2,8 . Though low nitrogen values could be explained by a heavier dependence on nitrogen poor browse material in bison diet 17 . Substantial incorporation of browse is reported by other Late Pleistocene bison paleoecology studies in North America based on stable isotopes and bison dentition wearing patterns 33,52 . Additionally, more canopy cover from prominent evergreens would lower the abundance of nitrogen in plant material consumed 29 . However, we cannot ignore that this could also be an effect of small sample size from only bison dentin for this temporal period.
Evergreen forests south of the Laurentide ice sheet were rapidly succeeded by other vegetative communities in the transition between the Pleistocene and the Holocene 49,53,54 . By the Early Holocene, new deciduous forest south of the evergreens formed and bordered along riparian areas while grasslands spread throughout the open landscape 54 . As the Great Plains quickly became dominated by prairie 49,54 bison migrated northward into the developing terrain and became plentiful at this time 53 . Yet only one data point from our sample assemblage falls within the Early Holocene boundaries. This single bison indicates an increase in C 4 vegetation incorporated into diet (Table 2), consistent with rising temperatures and the presence of more C 4 Chloridoideae grasses 51,55 . Additional evidence of an increase in C 4 grasses are recorded in bison living during the Early Holocene within present day Yellowstone National Park and the state of Nebraska 40,56 . Nitrogen levels are higher than in Late Pleistocene bison despite an increase in effective moisture in the Early Holocene 55 , further supporting the idea that bison selected more browse material from low light areas in evergreen forests during the Late Pleistocene.
The Middle Holocene climate is summarized as highly variable with an overall shift towards warmer, drier conditions and patchiness of resources 1 . Herbaceous plants fluctuated throughout the Middle Holocene, alternating between Poaceae and Ambrosia (ragweed) communities, indicating frequent changes in precipitation 50,57 . Bison from this sub-epoch exhibit the highest mean values for carbon and nitrogen as well as large variation in bone samples (Table 2, Figs 3, 4). Although the two bone samples from the Middle Holocene (Bison 22 and 98) coincide in both time and space, they exhibit large isotopic differences (Table 1, Figs 1, 3). This could be representative of the vastly fluctuating climate or a remnant of different feeding strategies among bison sexes 37 . In any case, we observe the greatest amount of C 4 vegetation in bison diet at this time for both tissue types (24%), indicating a trend towards longer and warmer growing seasons 39 . This is corroborated by bison isotopic values from the Eastern Great Plains during the Middle Holocene 11 .
Climatic conditions in the Northern Great Plains during the Late Holocene generally followed a cooling trend with increasing moisture up to modern day 26,55 . However, climate proxies provide evidence that severe arid conditions occurred at intervals throughout this time period 26 www.nature.com/scientificreports www.nature.com/scientificreports/ in bison samples from the Late Holocene but changes in vegetation type may lower the amount of nitrogen available in soils and dampen the signal of physiological stress 58 . A wide range of δ 13 C values and more C 4 plant material is recorded in Late Holocene bison (Table 2, Figs 3, 4), suggesting diverse vegetation utilized and a continued shift towards extended growing seasons 39,41 . Pollen records indicate that Poaceae increased in abundance during the Late Holocene and the first appearance of Salaginella densa (spikemoss) is documented in the Northern Great Plains 50 . The ground cover provided by spikemoss and its ability to persist in dry conditions provides protection from erosion and forage to subsist on during lean winter months. Its expansion likely increased foraging capacity and contributed to the immense presence of bison on the prairie during the Late Holocene.
Coinciding with the environmental changes that took place during the Late Holocene is the rise of a more complex human ecosystem throughout North America. How humans influenced landscapes, bison behavior, and available food supplies during the Holocene are currently not well understood but research suggests they had an active role in the Great Plains ecosystem 59 . Bison's past response to changing composition of habitat remains unclear but their ability to adapt and exploit a variety of resources is attributed to the species' long-term survival over other megafauna.
Finally, we compared isotopic values of modern bison dentin from TRNP with ancient bison from the Late Pleistocene through the Late Holocene. We observe comparatively low variability in modern bison stable isotopes (Table 2, Figs 3, 4). TRNP bison have a more restricted rangeland than their ancient counterparts as well as a presumably shorter window for sample collection. Modern bison are depleted in 13 C, similar to Late Pleistocene bison (Figs 3, 4), indicating a diet of 100% C 3 plant material (Table 2) despite different climatic conditions experienced by the temporal groups. Nitrogen values are notably lower than in other bison groups (Table 2, Fig. 4),  www.nature.com/scientificreports www.nature.com/scientificreports/ indicating no evidence of nutritional stress and adequate available moisture 27 . While low variation may be attributed to small sample size, other studies have found similar results in several living herds 25,39 . Tieszen (1994) 13 showed that the Wind Cave National Park bison herd in South Dakota had a diet with more C 4 plants but bison also contained a small amount of δ 13 C variability within the herd. The Catalina Island bison population also exhibits a comparable δ 13 C average and low variability 22 . Modern Yellowstone National Park bison exhibit low variability in δ 13 C values and a similar mean to the TRNP bison in this study despite their ability to cover much larger areas and complete substantial elevational migrations 37,60 . Whether this trend in low variability in modern bison diet is due to restricted rangelands and herd management practices or if it is a result of a narrowing in plasticity from the recent genetic bottleneck in bison history is still not clear. We would expect that if it were only due to the habitat restrictions imposed upon modern bison, we would observe more variability in Yellowstone herds.
Understanding the predecessors of present-day bison may unlock new views for reintroducing them more broadly to the North American landscape. These techniques are already influencing management decisions for European bison (Bison bonasus). European bison had to overcome similar environmental challenges as North American species during the terminal Pleistocene and are also predominately constricted to limited rangelands today 61 . Most habitat of modern European bison is forested but their morphological adaptations suggest they evolved in open grasslands 62 and then moved into woodland areas as the forests expanded and pressure from humans increased 29 . Isotopic studies of the ancient Eurasian steppe bison (Bison priscus) are informing conservation strategies for their ecological successor, Bison bonasus 63 . Several studies have found that Bison priscus did rely heavily on grazing, with more browse incorporated over time, as woody vegetation became more accessible 52,61 . This information makes the introduction of European bison to more open grassland habitats a plausible strategy for large scale restoration and is an example of the value of conservation paleobiology for current species management.
The North American landscape has been transformed dramatically during the last 250 years, and with few exceptions, bison are no longer allowed to migrate or range widely in localities where they currently exist. Further, the extreme population bottleneck experienced by bison at the end of the 19 th century has left the species with only a microcosm of the genetic toolkit that it once wielded for adaptation. Thus, both the resiliency of the species and the landscape it once inhabited have been altered in a manner unprecedented since the last ice age. We may expect that genetically isolated and spatially confined herds will be the most challenged by environmental fluctuations 64 . Range expansion efforts such as opening of new state, federal, and tribal lands to bison and establishment of conservation herds on private lands are already underway, and bison range is currently expanding, but only at incrementally small amounts in comparison to the native range of the species.
Despite limitations imposed on present day herds, isotopic relationships identified here have provided a unique glimpse at paleoecology that is relevant to current management of the species. Managers would benefit from a spatially and temporally expanded study of bison isotopic profiles throughout the known historic range of the species. An evaluation of historic genetic diversity of the species may also reveal pre-extirpation global and local herd level values that can be used in parallel with isotopic data to illuminate ancient bison ecology and inform management. The common occurrence of bison remains in archaeological collections curated in facilities throughout North America makes this effort feasible and demonstrates the value of such collections to present day management of species and systems.

Materials and Methods
Sample assemblage. Ancient bison bone and tooth samples were collected from four North Dakota museum collections: North Dakota Heritage Center and State Museum (State Historical Society of North Dakota), Knife River Indian Villages National Historic Site, University of North Dakota's Department of Anthropology, and University of North Dakota's Biology Museum. Samples included in this project were collected as part of twenty-one previously excavated archaeological sites in North Dakota and one previously excavated site in northern South Dakota (Fig. 1). The context of the archaeological sites encompasses a large temporal scale from the Late www.nature.com/scientificreports www.nature.com/scientificreports/ Pleistocene throughout the Holocene. Holocene sites are associated with each North Dakota Native American cultural tradition as defined by long standing archaeological schema including Plains Village (AD 1200-1800's), Plains Woodland (400 BC to AD 1700's), Plains Archaic (5,500 BC to 400 BC) and Paleo-Indian (11,500 to 5,500 BC). To avoid repeat sampling of bison individuals, we selected the right 3rd molar whenever possible. In other cases, we selected the molars or large premolars from the right side of the jawbone. Specimens were also chosen from different stratigraphic layers in the archaeological site context. While every effort was made to not repeat samples in the assemblage, there is a minute chance that repeats were made in some cases.
The integrity of skeletal elements was observed under stereo microscopy, revealing well preserved tissues at surfaces and on cross section of bones and teeth. Preservation quality was determined by a clear delineation between cortical and spongy material with little discoloration, compact tooth dentin and opalescent enamel ( Supplementary Fig. S3).
Modern bison teeth were obtained from animals (n = 5) culled for management purposes or those dying of natural causes within the North and South Units of TRNP. Bison range at TRNP includes two geographically separate units encompassing in total 28,542 ha the badlands of Western North Dakota. The Little Missouri River traverses park lands from south to north, and the landscape is characterized by stratified clay buttes capped with clinker and interspersed with lignite and fossils dating to the Paleocene Epoch. Annual precipitation is 38.1 cm, and vegetative communities include mixed grass prairie and sage (Salvia and Artemisia spp.) and other woody shrubs in uplands, cottonwood (Populus deltoides) galleries along riparian corridors, ash (Fraxinus pennsylvanica) groves in draws, and juniper (Juniperus spp.) stands along north facing slopes. Temperature varies widely, with means ranging from highs of approximately +29 °C in summer to lows of −18 °C in winter, with extremes sometimes ranging +43 to −43 °C.
At TRNP, bison forage alongside feral horses (Equus caballus), longhorn cattle (Bos taurus), elk (Cervis elaphus), pronghorn antelope (Antilocapra americana), mule deer (Odocoileus hemionus), white-tailed deer (Odocoileus virginianus), moose (Alces alces), sheep (Ovis canadensis), black-tailed prairie dogs (Cynomys ludovicianus), and other small mammals. Park lands are fenced, and in the absence of predators, bison are managed to prevent overgrazing through roundups, after which excess animals are transferred primarily to tribes. Herds have typically been allowed to range between 100-300 individuals in the North Unit and between 300-500 in the South Unit, in alignment with a forage allocation model and perceived social carrying capacity 48,65,66 . Research was approved by the University of North Dakota (UND) Institutional Animal Care and Use Committee (protocol number 1511-5) and the UND Institutional Biosafety Committee (registration number IBC-201511-007). All methods and experiments were performed within the guidelines and regulations for the use of experimental animal specimens.
All archaeological samples were selected from human derived contexts. This provides potential to access both bison ecology and human selection processes that resulted in the culling and ultimate deposition of individual skeletal elements in the archaeological record. We acknowledge the human induced bias in our sample. Taking this bias into account allows us to access a record of bison population dynamics otherwise unavailable as a result of the late 19 th century population bottleneck.
Sample preparation and isotope analysis. Sections of cortical bone or tooth dentin weighing 1-3 grams were cut from bison specimens using a band saw. Tooth dentin samples were taken in 0.5 to 1 inch contiguous pieces from dentin underlying enamel in the orientation of crown to cusp. Some tooth dentin portions included root material. Samples were then sent to the University of California Irvine (UCI) Keck Carbon Cycle Accelerator Mass Spectrometry Laboratory where they were decalcified in 1 N HCl and then gelatinized at 60 °C, pH of 2, and ultrafiltered to select for a high molecular weight fraction (>30 kDa). Aliquots of ultrafiltered collagen were measured on a Fisons NA1500NC elemental analyzer/Finnigan Delta Plus isotope ratio mass spectrometer to obtain δ 13 C and δ 15 N values at a precision of <0.1‰ and <0.2‰, respectively. Stable isotope measurements were determined on amino acid hydrolysate samples. Samples with C:N atomic ratios between 2.9 and 3.6 are indicative of well-preserved collagen 46,67 and were then measured for 14 C dating at the UCI Keck AMS facility using the methods outlined in Lohse et al. 11 . The error range for 14 C ages (BP) in this study is ±15-35 years. All 13 C to 12 C ratios were reported relative to the Vienna PeeDee Belemnite (VPDB) standard and all 15 N to 14 N ratios are reported relative to the Ambient Inhalable Reservoir (AIR) standard. 14 c calibration and temporal episodes. 14 C ages were calibrated using OxCal 4.3 68 and the IntCal13 curve 69 for the Northern Hemisphere. All bison sample ages are reported in calibrated years before present (cal BP) and are within a 95% confidence interval. The median of the confidence interval was used as the sample date for separation into temporal episodes. All δ 13 C values dated before 1800 AD were adjusted by -1.5‰ to account for the reduction in atmospheric CO 2 due to the increased burning of fossil fuels after the Industrial Revolution 13 .
The temporal range was split into episodes to allow comparisons between bison from the Late Pleistocene to modern. The episodes follow formal Holocene subdivisions recognized by the International Union of Geological Sciences (IUGS) and is based on data from Greenland ice cores, pollen records, lake sediments, and Global Stratotype Section and Points (GSSPs) 47 . The Pleistocene-Holocene boundary was marked at 11.7 thousand years ago. The Holocene was split into four episodes, Early Holocene (11.7 to 8.2 thousand years ago), Middle Holocene (8.2 to 4.2 thousand years ago), and Late Holocene (4.2 to 100 years ago) 47 . Modern bison are also considered an episode. δ 13 c and δ 15 n data analysis. We used generalized additive models (GAMs) to illustrate bison carbon and nitrogen isotopes over time with R statistical software and the package "mgcv" 70 . This method allows the estimation of isotope values between data points, providing a continuous view of bison isotopic fluctuations throughout the Holocene with a 95% confidence interval. GAMs were modeled with a gaussian distribution, an identity link function, and a smoothing parameter on time (k = 5).
The δ 13 C means of temporal episodes were also used to calculate the percentage of C 3 and C 4 grasses in bison using equation 1, modified from Carlson et al. (2018) 32  where δ 13 C collagen includes a 6.3‰ adjustment for trophic level fractionation specific to bison 16,66 , δ 13 C C4 = −12.5‰ and δ 13 C C3 = −26.5‰.

Data availability
All data generated for this study are included within this paper.