Abstract
Biological conservation depends on understanding and disentangling the effects of decadal- to centennial-scale dynamics from the millennial-scale dynamics documented in the fossil record. The American Southwest is expected to become increasingly arid over the next few decades and will continue to experience large-scale human land use in various forms. A primary question is whether the ecological fluctuations recorded over the past few decades fall outside the range of variation expected in the absence of recent land use and management. I excavated and quantified mammal diversity change in two fossil-bearing alcoves, East Canyon Rims 2 and Rone Bailey Alcove, located in San Juan County, Utah. AMS radiocarbon dates on 33 bone samples from these sites span ~4.4–0.5 kyr and shed light on pre-industrial faunal dynamics in the region over the course of environmental change. Localities with comparable small mammal diversity have not been reported from this region previously, so these deposits provide novel insight into Holocene mammal diversity in southeastern Utah. Taxa recorded in these sites include leporids, sciurids, perognathines, arvicolines, Onychomys, Cynomys, Dipodomys, Peromyscus, Neotoma, and Thomomys. Using temporal cross-correlation, I tested for a relationship between regional temperature and species richness, evenness, relative abundance, and rank abundance. I also tested for changes in the overall taxon abundance distribution and visualized faunal relationships among time bins using non-metric multidimensional scaling of relative abundance data. None of the measures of diversity tested here were correlated with temperature change through time except for relative abundance of leporids. Overall, these results suggest that climatic fluctuations of the magnitude preserved in these deposits did not significantly alter the small mammal community, nor is there evidence that the presence, then exodus, of Native Americans from the region significantly affected small mammals.
Keywords
Introduction
The biodiversity of the Colorado Plateau (CP; southwestern United States) is threatened by both land use and climate change impacts. This region is 42% federally owned land (less than 14% is privately owned), and over 80% of that land is multi-use, including mining, grazing, energy development (fossil fuel and solar/wind), off-road vehicle travel, timber harvesting, hiking, fishing, hunting, and other activities (Ernst and Prior-Magee, 2008). Understanding ecological dynamics on these multi-use lands is crucial for preservation of the ecosystem services and intangible benefits they provide (Fu et al., 2015; Lawler et al., 2014; Nash, 1967). While much information is becoming available on ecological dynamics that operate over decadal and shorter time scales, we still know little about the over-lying long-term dynamics that are important in conservation efforts (Conservation Paleobiology Workshop, 2012; Dietl and Flessa, 2011; Dietl et al., 2015; Hadly and Barnosky, 2009; Kidwell, 2015; Rick and Lockwood, 2013; Swetnam et al., 1999; Willis and Birks, 2006). A primary question is whether the ecological fluctuations recorded over the past few decades fall outside the range of variation expected in the absence of recent land use and management. For example, to what extent do current relative abundances, distribution, and associations of species reflect recent adjustments of species due to anthropogenic pressures versus natural fluctuations that typify ecological dynamics that play out over millennia? Such questions can only be answered by tracing ecological dynamics through thousands of years, using the natural experiments preserved in the fossil record (Conservation Paleobiology Workshop, 2012; Hadly and Barnosky, 2009; Kidwell, 2015).
Biological conservation on the CP depends on understanding and disentangling the effects of decadal- to centennial-scale dynamics, such as grazing and other human impacts, from the millennial-scale dynamics documented in the fossil record. The American Southwest is expected to become increasingly arid over the next few decades: conditions analogous to previous multi-year droughts, including the Dust Bowl, are expected to become the norm (Seager et al., 2007). Under the A2 (highest) greenhouse gas emissions scenario, precipitation is estimated to decline by around 66% by 2090 and temperatures are expected to increase ~1.5–3°C by 2041–2070 and ~3–5°C by 2070–2090 (Garfin et al., 2014).
The effects of climate on biota can take many forms: mammals can shift their geographic range, population abundance, physiology, body size, phenology, and so on, or, if a species’ environmental tolerances are greater than the amount of environmental change expected, they may remain observably unaffected. Assessing community-level fluctuations is one way to synthesize these many responses and compare the whole mammal community from one time period to another. Detailed paleoecological records from Quaternary deposits have been remarkably useful in characterizing these millennial-scale ecological dynamics in other regions (Anderson, 1993; Anderson et al., 1999, 2000; Barnosky et al., 2004; Betancourt, 1984; Blois et al., 2010; Graham and Grimm, 1990; Grayson, 2011; Hadly, 1996). Past studies have detected a strong relationship between faunal composition and climate: in Samwell Cave (California), Blois et al. (2010) show a correlation between warming climate and declining species evenness and richness from ~17 to 1.5 kyr; Terry et al. (2011) found that cave deposits in the Great Basin record proportional increases in ‘southern affinity taxa’ when climate warmed and dried over the last ~7.5 kyr; and the Porcupine Cave (Colorado) pit sequence spanning 1 ma to ~600 kyr revealed little community change due to climate change except that small herbivores were less diverse during warm interglacial periods (Barnosky et al., 2004). However, while small mammals are important indicators of climate and environment, they do not always track vegetation (Graham and Grimm, 1990) or temperature and precipitation predictably. For example, Lamar Cave in Yellowstone records the presence of mesic-adapted taxa during time intervals that are considered warm/arid (Hadly, 1996). At Mescal Cave in the northern Mojave Desert, Stegner (2015) found that although both Neotoma cinerea and Marmota flaviventris are considered boreal-adapted, N. cinerea disappeared at the end of the Pleistocene while M. flaviventris was present at the site during several thousand years of Holocene warming and aridification.
Few Quaternary fossil records that sample the small mammal community of the CP have been studied to date (Carrasco et al., 2005; Emslie, 1986; FAUNMAP Working Group, 1994; Mead, 1981; Mead and Phillips, 1981; Tweet et al., 2012), particularly in southeastern Utah, where only a handful of Quaternary vertebrate localities have been published. Most of these contain fewer than five taxa – primarily large-bodied species – and few specimens (Carrasco et al., 2005; FAUNMAP Working Group, 1994, 1996). Archaeological and fossil plant records from the Quaternary CP, in contrast, have been more-extensively studied and provide an important context for new vertebrate fossil data presented here. Holocene rockshelter deposits are common across the CP (Mead et al., 2003; Tweet et al., 2012) and contain abundant small mammal remains that can be used to track the communities through long time spans. To this end, I excavated and quantified mammal diversity change in two fossil-bearing alcoves located in San Juan County, Utah. These localities, East Canyon Rims 2 (ECR2) and Rone Bailey Alcove (RBA), contribute to our understanding of natural variation in this system by providing faunal data for a period of recent climate change – a transition from generally cool-wet to warm-dry that occurred within the last 5 kyr. Accelerator Mass Spectrometry (AMS) radiocarbon dates on 33 bone samples from these sites span ~4400 cal. BP to present and shed light on pre-industrial faunal dynamics in the region over the course of environmental change, most notably aridification. Localities with comparable mammal diversity have not been reported from this region previously, so these sites provide novel insight into Holocene mammal diversity in southeastern Utah.
Study region
The CP is a physiogeographic province in North America centered on the Four Corners Region where Utah, Colorado, New Mexico, and Arizona meet. It is an arid region with high topographic relief – from ~900 to 4300 m – pocked with isolated laccolithic mountains and creased with deep canyons formed by the Colorado River and its tributaries (Barnes, 1993; Chronic, 1990). The CP is flanked to the west by the Great Basin, to the east by the Rocky Mountains, and to the south by the Sonoran and Chihuahuan Deserts. The western border is formed by a series of mountain ranges – the Uinta Mountains, Wasatch Range, and Fishlake and Aquarius Plateaus (Rowe, 2007) – while the Mogollon Rim in Arizona and the Rio Grande Rift Valley in New Mexico form the southwestern and southeastern borders, respectively (Foos, 1999).
Localities
RBA and ECR2 are large, horizontally shallow alcoves in the Slickrock Entrada Sandstone cliff face of Rone Bailey Mesa, Canyon Rims Recreation Area (BLM), San Juan County, Utah (Figure 1a–c). ECR2 is 1830 m (±5 m) in elevation, southeast-facing, and roughly ~50 m across, ~30 m high at the mouth, and ~20 m deep (Figure 1c). RBA is 1905 m (±5 m) in elevation and faces roughly southwest. RBA (~10 m high by ~3 m wide by ~5 m deep) is contained within a much larger alcove, and an apron of sediment roughly 10–15 m high extends from a maximum height at the midden, spilling southwest to the floor of the larger alcove (Figure 1b). Quaternary eolian sediments have accumulated in both ECR2 and RBA and buried bone and plant material collected by roosting raptors (as evidenced by raptor pellets on the surface of the sediments) and woodrats (genus Neotoma). Woodrats remain active at both sites today, and mammalian carnivores also may have played a role in bone accumulation in these middens.
(a) Map of localities and geographic landmarks mentioned in the text; regional map is inset, the dashed line indicates the rough outline of the Colorado Plateau. (b) RBA from approximately west; arrow indicates location of woodrat midden deposit. (c) ECR2 from approximately south; arrow indicates location of woodrat midden deposit.
The vegetation surrounding these sites today is a patchy mixture of sagebrush-dominated areas, open perennial grassland, and pinyon–juniper woodland. Cooler and more mesic pockets harbor gambel oak (Quercus gambelii), barberry (Mahonia fremontii), squawbush (Rhus trilobata), Utah service berry (Amelanchier utahensis), and spruce (Picea englemannii). At ECR2, invasive Russian thistle (tumble weed, Salsola spp.) dominates the middle of the drainage where a primitive road passes through and forms dense patches up to approximately a meter high in places. Russian thistle is present near RBA as well, but it is sparser and does not produce the monotypic thickets seen near ECR2. The area around Rone Bailey Mesa has been grazed by cattle and horses for over a century (Lavender, 1964).
Climate and vegetation history of CP
Today, the CP marks the geographic transition between summer-wet (summer monsoon) to the south and summer-dry to the northwest (Anderson et al., 2000). Southeastern Utah currently experiences summer monsoons, but the monsoon boundary shifts through time in response to temperature, snow pack in the Rocky Mountains, location of the Jetstream, and other factors (Anderson et al., 2000). Because the CP has high topographic variability – deep canyons and high peaks – typical precipitation and temperature are spatially heterogeneous. However, most of the region receives between 4 and 12 inches of precipitation annually and experiences minimum January temperatures around −10°C to −4°C (13–25°F) and maximum July temperatures between roughly 29°C and 38°C (84–100°F) (30-year normals from 1981 to 2010) (PRISM Climate Group, 2015).
The CP was cooler and more mesic during the last glacial period of the Pleistocene, before ~14 kyr, at which time it was characterized by juniper woodland and sagebrush between ~600 and 1500 m; pinyon, limber pine, and Douglass fir between 1500 and 1800 m; Engelmann spruce/juniper forest around 2100 m; and spruce/pine forest above 2700 m (Anderson et al., 2000; Cole, 1990). Since the last glacial maximum, the dominant tree and shrub species – including Abies concolor (white fir), Artemisia spp. (sagebrush), Atriplex confertifolia (shadscale saltbush), Cercocarpus intricatus (mountain mahogany), Ephedra spp. (Mormon tea), Juniperus spp. (Juniper), Opuntia spp. (prickly pear), Picea spp. (Spruce), Pinus ponderosa (Ponderosa pine), Yucca angustissima (narrow-leaf yucca), and many others – have generally moved 700–900 m higher in elevation and 400–700 km up river (Cole, 1990). At the end of the Pleistocene, middens from the Abajos and Comb Ridge (respectively ~45 and ~85 km south of Rone Bailey Mesa) record a mixture of xeric- and mesic-adapted plants; the mesic components (e.g. subalpine conifers) of the flora disappeared from lower elevations, but modern dominant trees – pinyon (Pinus edulus and Pinus monophylla) and ponderosa (P. ponderosa) – did not appear at those elevations until the mid-Holocene (~7–3 kyr) (Betancourt, 1984; Coats et al., 2008). Pinyon has been expanding from south to north across the CP since the end of the glacial period (Coats et al., 2008).
With increasing elevation in this arid region, precipitation increases and temperature decreases on average (Betancourt, 1984). The fossil localities are between ~1830 and 1900 m in elevation, so the plant community likely progressed from pinyon, Douglas fir, and other conifers at the end of the Pleistocene, to juniper woodland and sagebrush with smaller pockets of ponderosa and spruce over the course of the Holocene, to an assemblage dominated by sagebrush and grasslands with large pockets of pinyon–juniper today. Preliminary analysis of pollen from RBA has identified spruce in the deepest excavation level, around 4400 cal. BP, and juniper throughout the deposit.
Climatically, the early Holocene was cooler than today, but more mesic than the late Glacial because the summer monsoon was strengthened (Weng and Jackson, 1999); this is also when the modern monsoon boundary was established (Betancourt, 1984). This cool mesic period gave way to an arid and warm mid-Holocene, from about 8.5 to 6 kyr (Reheis et al., 2005; Weng and Jackson, 1999). From ~6 to 3 kyr, cool-wet conditions returned (Reheis et al., 2005) and fossil plant evidence from the Abajos suggests higher effective moisture before 3 kyr (Betancourt, 1984). At the end of this period, around 3120 cal. BP, Maize agriculture arrived in this area as indicated by pollen in the Abajos records (Betancourt and Davis, 1984). Analysis of eolian and alluvial deposition in Canyonlands National Park suggests that from 2 kyr to the present drier conditions set in, as evidenced by greater mobility of dune sand (Reheis et al., 2005). In contrast, flood deposits from Utah and Arizona suggest drier periods from 3.6 to 2.2 kyr and 0.8 to 0.6 kyr, contrasted with peak flood magnitude and frequency from 1.1 to 0.9 kyr and again after 0.5 kyr (Ely, 1997). These peak floods are linked to increases in winter storms and dissipating tropical cyclones rather than summer precipitation or monsoon (Ely, 1997). Tree ring data record a mesic period from ~900 to 820 BP, followed by a series of at least six muti-decadal (>20 years) droughts – so-called megadroughts – that took place from ~1100 to 350 BP (Benson and Berry, 2009). The first of these megadroughts occurred between 1087 and 1066 BP (Benson and Berry, 2009). Human settlements in the northern San Juan Basin (southeastern Utah and southwestern Colorado) were largely abandoned by around 650 BP (Benson and Berry, 2009). Whether this exodus is attributable directly (e.g. through crop failure) or indirectly (e.g. through social collapse and warfare possibly due to food shortage) to climatic changes is still debated, but climate-induced stress is a predominant theory today (Benson and Berry, 2009; Stahle and Dean, 2010).
In the last two centuries, large-scale livestock grazing on the CP has transformed aspects of this region – among many impacts, livestock break up soil crusts and destabilize soils, leading to increased dust, soil loss over time, and exotic plant species invasion (Belnap, 2003; Belnap et al., 2009; Miller et al., 2011); grazing also changes the plant structure including dominance by weedy invasive species (Belnap et al., 2009), increases wild fire potential (Belnap et al., 2009), reduces carbon storage and vegetation cover (Fernández et al., 2008; Jones, 2000; Miller et al., 2011), and leads to destruction of riparian habitats and waterways as livestock congregate around water sources (Eldridge and Whitford, 2009). While the midden data addressed here do not provide clear records for faunal change between 250 BP and present, when livestock grazing began and then intensified, these localities elucidate what we might expect from the small mammal community in response to climate if grazing has no effect.
Methods
Excavation and identification of vertebrate fossils
Middens like RBA and ECR2 provide high-fidelity records of diversity through time. These sites generally sample the small vertebrate community in ~7 km radius around the site (Feranec et al., 2007; Hadly, 1999; Porder et al., 2003). Because ECR2 and RBA are ~4.7 km apart, thereby sampling an overlapping area, and they have similar taphonomic vectors, they can be treated as representative of the same community and the stratigraphic levels from the two sites can be combined to construct a chronological sequence that extends from ~4.4 to 0.5 kyr. Taphonomic studies of woodrat- and owl-generated deposits reveal that they record relative abundance and diversity with high accuracy (less than 1% mismatch between modern middens and the living local community) (Terry, 2008, 2010a, 2010b). Deposits like these are also spatially integrated, sampling small mammals from all microhabitats surrounding the deposits (Terry, 2010a). Time-averaging on the scale of decades to centuries is actually an asset to studies of diversity change through time because time-averaged middens sample rare species which are often missed by other survey techniques (Barnosky et al., 2004; Hadly, 1999; Terry, 2010b), smooth the effects of annual or decadal population booms and busts which are common in small mammals (Dickman et al., 2010; Terry, 2008; Whitford, 1976), and these sites accurately reflect the magnitude and direction of faunal baseline shifts (Terry, 2010b).
RBA and ECR2 were identified as promising deposits in 2012 and 1 m2 test pits were excavated in each site in 2013. Excavation levels were an arbitrary depth based on sediment volume (ranging from 2.3 to 13.5 cm in thickness; Supplemental 1, available online) because sediments were extremely unconsolidated sand mingled with roof fall (large sandstone blocks), and there was no apparent natural stratigraphy. Sediments were screened (1/2, 1/4, and 1/16 inch2) and bone was picked from the 1/2- and 1/4-inch2 screens in the field. All material collected on the 1/16-inch2 screen was bagged, taken to the University of California Museum of Paleontology (UCMP), and processed for this study.
Fossil material was separated from matrix by hand or with forceps. Material was initially sorted into rodent scats, insects, plant macrofossils, and vertebrate bone. Using morphological criteria, bone was identified to class (Aves, Reptilia, Amphibia, and Mammalia), and mammalian bone was identified to species when possible (primarily craniodental remains), or to family or genus (in the case of postcranial or otherwise non-diagnostic fossils), by direct comparison with modern skeletal specimens of known species held at the Museum of Vertebrate Zoology (MVZ) and UCMP or using published descriptions from primary paleontological and mammalogical literature. All specimens are curated at the UCMP. See Supplemental 2 (available online) for Systematic Paleontology.
Radiocarbon dating
I radiocarbon dated 33 bone collagen specimens (16 from ECR2 and 17 from RBA; Table 1) using AMS dating techniques at the Center for Accelerator Mass Spectrometry, Lawrence Livermore National Laboratory (Livermore, California, US). Seven of these specimens were subsampled and dated twice to confirm that results were consistent. Bone specimen preparation follows the modified Longin method described by Brown et al. (1988) for collagen extraction and follows methods outlined in Vogel et al. (1987) for converting CO2 into graphite for AMS analysis. AMS 14C results include a matrix-specific background correction and an estimate of the δ13C value of the material (Stuiver and Polach, 1977). To convert radiocarbon years BP to calendar BP, I used the IntCal 13 calibration curve (Reimer et al., 2013) in OxCal Online, version 4.2 (Bronk Ramsey, 2009).
AMS radiocarbon and calibrated ages of bone collagen samples from ECR2 and RBA.
AMS: Accelerator Mass Spectrometry; CAMS: Center for Accelerator Mass Spectrometry; UCMP: University of California Museum of Paleontology; RBA: Rone Bailey Alcove.
Specimens that were subsampled and dated twice.
Faunal data for excavation levels with overlapping radiocarbon age ranges (Figure 2) were combined, and I developed a chronology of six time bins with different degrees of time averaging, dating from 4410 to 501 cal. BP.
Calibrated age (cal. yr BP) of bone collagen samples from ECR2 and RBA. Diamonds indicate mean cal. yr BP, and vertical lines indicate the age range. Brackets indicate the excavation levels that were combined to produce time bins (letters).
Sample standardization
Taxon abundance was quantified as number of individual specimens (NISP) of craniodental material identified to family (Leporidae, Sciuridae), subfamily (Perognathinae, Arvicolinae), or genus (Onychomys, Cynomys, Dipodomys, Peromyscus, Neotoma, and Thomomys). As compared with minimum number of individuals (MNI), NISP has been demonstrated to be a less-biased indicator of relative importance of a taxon in the community (Blois et al., 2010; Grayson, 1978). Species-level identifications within these taxa are often only possible by comparison of skeletal elements that are seldom well represented in the fossil record. Analysis at taxonomic levels above the species level allows incorporation of more of the fossil data represented in these deposits, improving abundance estimations by increasing sample size. Furthermore, analysis at the genus, family, and subfamily levels provides the requisite information needed to recognize the relative changes in the small mammal taxa that fill different functional roles in the ecosystem (Barnosky, 2004; Barnosky et al., 2004; Blois et al., 2010; Grayson, 2000; Terry et al., 2011) and which can be used as appropriate metrics of the environmental state (Badgley and Fox, 2000; Barnosky and Shabel, 2005; Barnosky et al., 2004; Stegner and Holmes, 2013).
I calculated relative abundance as NISP of the taxon divided by the total number of specimens for each time bin. I then estimated 95% confidence intervals around each relative abundance estimate using Goodman’s (1965) simultaneous confidence interval calculations (Calede et al., 2011; McHorse et al., 2012), implemented in the DescTools package in R (Signorell et al., 2015).
Species richness is correlated with sample size, so I standardized richness using shareholder quorum subsampling (SQS) at a coverage level of 0.6 (Alroy, 2010). While the goal of traditional rarefaction is uniform sampling, the purpose of SQS is ‘fair’ sampling that is reflective of the true standing diversity. In SQS, sample size is not fixed; rather, SQS fixes coverage, the proportion of the entire frequency distribution represented by the species in the resample. Rarefaction dampens the richness of very diverse samples, so relative diversity estimates are not linear, whereas in SQS a site that is twice as diverse as another has an SQS richness value that is twice as large (Alroy, 2010). Estimated richness values produced by SQS appear low relative to raw richness, but are linearly proportional to one another.
Estimation of confidence intervals is one way to correct for differences in time bin duration, but longer time bins inevitably have more opportunity to sample more species, and low sample size, as in bins F and A, can lead to a randomly biased NISP and abundance distribution structure. I tested for a correlation between duration and both raw and standardized richness per time bin and found no relationship between the two (raw richness: Pearson’s product moment correlation = 0.017, p = 0.97; standardized richness: Pearson’s product moment correlation = –0.276, p = 0.60).
All analyses were performed using the R Project for Statistical Computing version 3.2.0 (R Core Team, 2015) and with the vegan package (Oksanen et al., 2015), DescTools package (Signorell et al., 2015), and SQS function (Alroy, 2010).
Analysis of diversity change
I assessed diversity change through time in several different ways. First, I used pairwise Fisher’s tests to compare the relative abundance distribution from one time bin to the next (with a Monte Carlo p-value simulation, 1000 replicates, and both with and without a Holm p-value correction). Second, I calculated Hurlbert’s probability of interspecific encounter (PIE) for each time bin and compared across time bins. PIE is a measure of species evenness and is independent of sample size, unlike many other commonly used diversity metrics, for example, Shannon and Simpson indices (Blois et al., 2010; Gotelli and Ellison, 2013; Hurlbert, 1971). Because PIE is not biased by sample size, resampling is not necessary for estimating the uncertainty in the metric. Instead, I followed the equation presented in Davis (2005) for computing variance in PIE. Because time bins are likely to be serially autocorrelated and therefore non-independent, I compared the overlap in variance to determine whether PIE was significantly different among time bins (Blois et al., 2010).
I visualized faunal relationships among time bins using non-metric multidimensional scaling (NMDS) of relative abundance data. Unlike other ordination techniques (e.g. principal components or principal coordinates analysis), NMDS does not preserve true distances between objects (time bins or taxa in this case) and instead uses rank order of the distances between objects to arrange them in multidimensional space. This technique places more different objects further apart and more similar objects closer together, which is not always the case in other ordination techniques (Gotelli and Ellison, 2013).
Rank-order abundance was scaled from 1 to 10, to account for time bins where not all species were represented or when some taxa had the same number of specimens. I made pairwise comparisons of rank-order abundance between time bins using Spearman rank correlation tests with and without a Holm p-value correction. The Spearman coefficient ranges from −1, indicating an inverse relationship, to 1, indicating perfect correspondence (Terry, 2010a). To determine whether particular taxa were consistently higher or lower in rank than the others, I compared the ranks of taxa using pairwise Mann–Whitney U-tests with and without Holm p-value corrections.
Climate correlation
To evaluate the role of climate across all six time bins, I used reconstructed mean July temperature (MJT) for the Southwestern United States from Viau et al. (2006), downloaded from the National Climatic Data Center (http://www.ncdc.noaa.gov). These reconstructions are based on 752 fossil and 4590 modern pollen records from across the continent and estimate MJT in 100-year intervals from 14 kyr to the present (Viau et al., 2006). Following the methods described in Terry and Rowe (2015) and Rowe and Terry (2014), I fit the climate data with a spline smoothing function (R function spline: stats, method: natural) and then analytically averaged the climate data over the period of time encompassed by each time bin. I then used cross-correlation to test for a relationship between climate and species richness, evenness, and relative and rank abundance of each taxon. The significance of the correlations was tested using permutation tests (500 iterations). To evaluate lagged correlations, I shifted the climate data by 10-year increments, up to 200 years, each time re-averaging the fitted climate data and recalculating the correlation (Terry and Rowe, 2015).
Results
Radiocarbon chronology
Radiocarbon chronologies for each site show that most dates were stratigraphically consistent, with levels ranging from 143 (ECR2 level 3) to 841 (RBA level 6) years of time averaging (Table 1, Figure 2). Dates presented in Figure 2 are derived from the most likely age for each specimen. However, there is a range of possible calibrated cal. BP for a given 14C date (Figure 2); I used the probability distribution associated with each date and the overlap in possible calibrated ages to determine which excavation levels could reasonably be considered non-overlapping in age. The three older time bins pertain to RBA, while the younger are from ECR2. Time bin A, the oldest, represents excavation level 6 from RBA and has peak age probabilities at 3699.5, 3639.5, and 4184.5 cal. BP. RBA level 5 was not dated, but dates for levels 6 and 4 are consistent and non-overlapping, so here I assume that stratigraphic mixing in these older layers is either not present or very minimal and that level 5 can be treated as a separate faunal unit (time bin B). However, additional radiocarbon dates are necessary to confirm this assumption. Level 4 encompasses time bin C with peak age probability at 2424.5 cal. BP. Time bin D represents ECR2 excavation levels 8–10, with peaks at 1024.5, 1169.5, 1179.5, 1194.5, 1214.5, 1289.5, 1529.5, and 1609.5 cal. BP. Time bin E represents ECR2 levels 4–6 and peaks at 794.5, 789.5, and 959.5 cal. BP. Time bin F is represented by ECR2 level 3. The oldest date from ECR2, 4970–4844 cal. BP, is from excavation level 3 and is ~2.6 kyr older than the next oldest date (from level 8). If this is excluded, bin F/level 3 has peaks at 524.5 and 544.5 cal. BP. While there is no evidence that this particularly old date is inaccurate, either as a result of contamination in the lab or in the deposit, these kinds of contamination are extremely difficult to diagnose. However, all other dates for ECR2 are consistently younger than 1704 cal. BP, and so this date of nearly 5 kyr in the middle of the deposit has a few important implications. First, it suggests that stratigraphic mixing is present to a certain degree in this site, but that it occurs rarely because dates from other levels tend to be consistent and arranged in a predictable chronological order. ECR2 level 1 and RBA levels 1 and 2 both have wider age ranges and show more mixing than the deeper/older excavation levels. Those levels have not been included in these analyses because the magnitudes of time averaging they encompass entirely overlap deeper excavation levels. ERC2 level 2 and RBA level 3 were not dated and were also excluded because it is not known whether they are highly mixed, and it appears that younger levels may be more mixed than older levels. Bioturbation by people, cattle, and rodents is one possible explanation for the presence of older bones in top excavation levels. Second, while there is some mixing in these sites, dates generally get older with increasing depth, suggesting it is reasonable to assume that the majority of the material in each level has not moved dramatically up or down in the deposit. Combining neighboring excavation levels into single analysis units helps to strengthen this assumption. However, interpretation of the age of individual specimens is not possible without a date on the specimen itself.
Relative and rank-order abundance
Abundance of taxa fluctuates across the record (Table 2, Figure 3), but 95% confidence intervals demonstrate that there is no significant difference in relative abundance through time within each taxon and few differences among taxa (Figure 3). In time bin A, observed relative abundances of arvicolines and Thomomys were at their highest, while leporids, perognathines, Cynomys, and Onychomys were not sampled at all. During B, leporids, perognathines, Cynomys, and Onychomys first appeared, and perognathines were at their most abundant across the record. Dipodomys also increased in abundance, while sciurids, arvicolines, Peromyscus, Neotoma, and Thomomys all declined. In bin C, leporids, Cynomys, Neotoma, Onychomys, Peromyscus, and arvicolines all increased, while sciurids, Thomomys, Dipodomys, and perognathines declined. Leporids and Cynomys continued to increase in bin D, along with slight increases in Neotoma, Onychomys, Peromyscus, and sciurids. Dipodomys and Thomomys continued to decline, Thomomys reaching its lowest abundance across the record, and arvicolines also decreased. Cynomys was again absent during bin E; Neotoma declined while Onychomys and Peromyscus slightly increased. After relatively low abundance in the preceding time bin, Thomomys sharply increased in E. Perognathines, Dipodomys, and sciurids also increased while leporids slightly declined. Bin F is characterized by extremes: leporids, Thomomys, and Neotoma were at their highest abundance, but Onychomys, Cynomys, and arvicolines were entirely absent. Dipodomys declined slightly, and Peromyscus, sciurids, and perognathines maintained low abundance.
Number of individual specimens (NISP) (above) with relative abundance (taxon NISP/total NISP) and 95% Goodman’s simultaneous confidence intervals (italics) (below); uncorrected species richness and standardized (SQS) richness.
SQS: shareholder quorum subsampling.
Relative abundance of each taxon through time. Vertical bars are Goodman’s simultaneous 95% confidence intervals. The order of time bins is illustrated for Dipodomys only; the same order applies for all taxa.
The rank-order abundances are all positively correlated across time bins (Figure 4, Table 3). Without the Holm p-value correction for multiple statistical tests, the rank-order abundances of B, E, and F are significantly correlated with one another; additionally, C is correlated with A and E, and D is correlated with F. All of these comparisons have correlation coefficients over 0.60. With a Holm correction, only E versus C and F versus E are significantly correlated. The Mann–Whitney U-tests (Table 4, Figure 4) show that, without a p-value correction, Dipodomys and Thomomys are significantly more rank abundant (more common taxa have lower rank abundance values; that is, rank of 1 is the first most common taxon) than all taxa except for leporidae and Neotoma (only Dipodomys is significantly higher in rank than Neotoma). Arvicolines are significantly less rank abundant than all other taxa except for Onychomys, perognathines, and sciurids; Cynomys is significantly less rank abundant than all taxa except Onychomys, perognathines, and arvicolines. No comparisons are significant when a Holm p-value correction is used.
Boxplot of rank-order abundance ranks across time bins for each taxon.
Correlation coefficients from pairwise Spearman rank correlation tests on rank-order abundance of taxa for each time bin (above); Fisher’s test uncorrected p-values, for pairwise comparison of abundance distribution among time bins (below, in italics).
Statistical significance without a Holm p-value correction.
Significance with a Holm correction.
The p-values from pairwise Mann–Whitney tests.
Significance without a Holm correction. No pairwise comparisons were significant with a Holm correction. Mean rank for each taxon is in parentheses.
Diversity dynamics
The abundance distribution of bin D was significantly different (Fisher’s test without p-value correction) from all time bins except E (Table 3). Bins B and F were also significantly different from one another. NMDS (Figure 5) illustrates the relative degree of compositional similarity among time bins, although no time bins form tight groups with one another. Sample size is low in bins A (n = 13) and F (n = 39), which could have a meaningful impact on the results of the Fisher’s tests and on the arrangement of time bins and taxa in the NMDS.
Plot of non-metric multidimensional scaling of relative abundance.
PIE (Figure 6) was relatively high in time bins A (4410–3569 cal. BP), C (estimated 2700–2350 yr. BP), E (1054–735 cal. BP), and F (644–501 cal. BP). It reached a peak during bin D (1704–959 cal. BP), but was significantly lower than other times during bin B (3569–2700 cal. BP). Variance is non-overlapping for all time bins, except A and F overlap.
PIE index for each time bin; vertical black bars and numbers are the variance associated with each PIE value. Letters indicate time bin and horizontal black bars are the age range (cal. yr BP). Time bin B was not dated, so the age range is unknown. The gray line and points are reconstructed mean July temperature anomaly from Viau et al. (2006).
Climate correlations
There were no significant correlations between MJT and diversity using temporal cross-correlation except for relative abundance of leporids, which had the highest correlation coefficient at lag = 0 (r = 0.83, p = 0.026). Because sample sizes are low, particularly in bins A and B, there is high degree of uncertainty around the relative abundance of each taxon, illustrated in Figure 3. Additional specimens could change these results by changing the relative abundance distribution for each time bin. Additional specimens are more likely to change the rank abundance for less common taxa (ranks above 3) because the magnitude of difference in NISP is smaller among the less common taxa than among the most common taxa. In short, larger sample size is unlikely to change the fact that Dipodomys, Thomomys, and leporids are generally the most common taxa.
Discussion
Overall, neither diversity nor abundance in ECR2 and RBA correlates with the climate record assessed here. There were no significant differences in abundance of taxa throughout the record, from ~4.4 to 0.5 kyr. Most of the taxa present in these deposits are tolerant or characteristic of xeric ecosystems, and only one taxon present in the fossil record is not found in the immediate vicinity of the caves today (voles, genus Lemmiscus and Microtus).
Species of Microtus are typically found in grasslands and meadows, often associated with montane and/or riparian ecosystems (Kays and Wilson, 2009). In the fossil record, Microtus tends to be associated with cooler and more mesic times (Barnosky, 2004; Hadly, 1996; Hadly and Barnosky, 2009). In contrast, the sagebrush vole, Lemmiscus curtatus, is found in more arid brushy environments – usually sagebrush or rabbit brush (Carroll and Genoways, 1980). L. curtatus is restricted to the Great Basin today but is found in late Pleistocene records from across the CP (Mead et al., 2003; Murray et al., 2005). Teeth of L. curtatus are present in ECR2 level 4 (time bin E, 1054–735 cal. BP) and RBA level 4 (time bin C, estimated 2700–2350 ybp). To my knowledge, these are the youngest records of L. curtatus on the CP (Bell and Glennon, 2003; Carrasco et al., 2005; FAUNMAP Working Group, 1994, 1996; Mead et al., 2003; Neotoma Paleoecological Database, 2013): other occurrences of this species on the CP in the last 40 kyr are Sheep Camp Shelter (late Glacial/Holocene; Gillespie, 1985), Screaming Neotoma Cave (29–25 kyr; Bell and Glennon, 2003; Mead et al., 2003), and Isleta Cave No. 2 (late Glacial/Holocene; Harris and Findley, 1964).
As MJT increased, leporid relative abundance also increased. However, the 95% confidence intervals around relative abundance of leporids overlap for all time bins, so it is possible that the observed increase is a sampling artifact. Sylvilagus and Lepus are both ubiquitous across the CP and in western North America generally and, as genera, are tolerant of a wide range of habitats (Kays and Wilson, 2009). In a study of Chihuahuan Desert Lepus californicus, Hernández et al. (2011) found a correlation between leporid abundance and both precipitation and plant productivity, and conclude that food supply has a strong effect on L. californicus abundance. In contrast, Lightfoot et al. (2010) also assessed leporid abundance dynamics in the Chihuahuan Desert but found no correlation between Sylvilagus audubonii or L. californicus abundance and precipitation or productivity. Rather, Lightfoot et al. (2010) observed higher diversity of both species in black gramma grassland as compared with creosote bush shrubland. L. californicus favors areas with sparse grass and vegetation cover, including heavily grazed areas, but avoids tall grasses and forests (Best, 1996). In contrast, the species of Sylvilgaus present in southeastern Utah today (S. audubonii and Sylvilagus nuttali) are more commonly associated with brushy areas (Chapman, 1975; Chapman and Willner, 1978). It is not generally possible to distinguish Lepus from Sylvilagus using morphological criteria or size of loose teeth (Robinson and Matthee, 2005; Stegner, 2015), so it is unclear whether the increase in leporids in RBA and ECR2 is driven by one genus or whether both became more abundant through time. Climate was almost certainly becoming more arid at these sites after about 2 kyr (Cook and Krusic, 2004; Pederson et al., 2011; Reheis et al., 2005) and local changes in the plant community were likely taking place, which may have favored leporids. Identification of the pollen from these sites would shed light on this question.
A second pattern is qualitatively evident in the relative abundance data (Table 2): when Thomomys have lower relative abundance, Neotoma and Peromyscus are higher in abundance. Also, as Dipodomys declines through time, leporids increase. However, Dipodomys and Thomomys remain the most common (first and second most abundant, respectively) throughout the record, except bin D (leporids are first in rank abundance, Dipodomys are second, and Thomomys are fifth) and bin F (Thomomys are first, leporids are second, and Dipodomys are third). Terry et al. (2011) found that granivores tended to increase in abundance as climate warmed: of the two granivores at ECR2 and RBA, Dipodomys actually decline in abundance as climate warmed, and perognathines appear to fluctuate at low abundance independently of climate. Abiotic affinities do not seem to be related to this relative abundance pattern: Dipodomys are physiologically adapted to dry environments (Vimtrup and Schmidt-Nielsen, 2005) and are found in a range of arid and semi-arid habitats (Garrison and Best, 1990; Kays and Wilson, 2009); Dipodomys ordii, the only species found in northern San Juan County today, is almost always associated with sandy soils (Garrison and Best, 1990). The two species of Thomomys present in Utah today are Thomomys talpoides, found generally in the mountains and high valleys of Utah (Durrant, 1952), and Thomomys bottae, generally found in agricultural fields (such as alfalfa), lower valleys, and desert mountains (Durrant, 1952; Verts and Carraway, 1999). However, both have extremely wide geographic, elevational, and ecological ranges (Durrant, 1952; Jones and Baxter, 2004; Verts and Carraway, 1999), and so their present-day distribution in Utah is not necessarily a good indicator of their fossil distribution or habitat preferences.
Peromyscus and Neotoma are also fairly generalist in their climatic requirements and are dietary generalists as well (Blois et al., 2010; Sorensen et al., 2004). Certain species within these genera are specialists, however – for example, Peromyscus truei (pinyon mouse) is found associated with pinyon (Hoffmeister, 1981), and Neotoma stephensi is associated with juniper or, occasionally, other conifers (Sorensen et al., 2004). Finer taxonomic resolution of the fossil data would certainly aid in interpretation of abundance dynamics, for example, by revealing intra-generic turnover, but isolated teeth of Peromyscus and Neotoma cannot be identified to species using standard morphological criteria (Grayson, 1988; Hooper, 1957; Mead and Phillips, 1981; Mead et al., 2003; Repenning, 2004; Stegner, 2015). Using the dental morphological criteria developed by Hooper (1957), I attempted to identify Peromyscus teeth from ECR2 and RBA to species. However, almost all teeth were morphologically aligned with at least two of four species present in this region today: Peromyscus maniculatus, Peromyscus crinitus, Peromyscus boylei, and P. truei (Supplemental 3, available online).
The few sciurid specimens from ECR2 and RBA are almost exclusively loose cheek teeth, and the high degree of wear on these teeth makes identification to genus challenging, with the exception of Cynomys which is more hypsodont than other sciurids and larger than all other sciurids except Marmota. Today, the sciurids present in southeastern Utah include the burrowing Cynomys gunnisoni, Ammospermophilus leucurus (white-tailed antelope squirrel), Otospermophilus variegatus (rock squirrel), and Xerospermophilus spilosoma (spotted ground squirrel); largely ground-dwelling Tamias mimimus (least chipmunk), Tamias quadrivitattus (Colorado chipmunk), and Tamias rufus (Hopi chipmunk); and tree-dwelling Sciurus aberti (Abert’s squirrel). All of these species, with the exception of S. aberti, are found in arid desert and scrubland habitats today (Belk and Smith, 1991; Best et al., 1994; Burt and Best, 1994; Kays and Wilson, 2009; Nash and Seaman, 1977; Oaks et al., 1987; Streubel and Fitzgerald, 1978; Verts and Carraway, 2001). I have trapped both Tamias and Ammospermophilus in patches of pinyon–juniper woodland near these sites. S. aberti is present in the Abajos today (Schaefer, 1991) but not outside montane, forested regions (Kays and Wilson, 2009) – the presence of S. aberti in the fossil record of ECR2 and RBA might indicate cool/mesic climatic conditions. The sciurid teeth recovered from ECR2 and RBA are likely too small to be Sciurus, but further research will be necessary to confirm. The Cynomys recovered from RBA and ECR2 are white-tailed prairie dogs, subgenus Leucocrossuromys, which still persist less than 1 km southeast of ECR2 today. Cynomys, even on a subgenus-level, is unlikely to be a good indicator of climate: prairie dogs can persist in semi-desert (C. (L.) leucurus), grasslands and prairies (C. (L.) parvidens, C. (C.) ludovicianus), as well as montane meadows (C. (L.) gunnisoni, C. (L.) leucurus) (Kays and Wilson, 2009).
Two additional arid-adapted taxa are present throughout the deposits in relatively low abundance. Today, perognathines – genus Perognathus and Chaetodipus – are generally found in arid or semi-arid grasslands, desert scrub, or shrub steppe, but several species also occupy woodlands and chaparral (Kays and Wilson, 2009). Only Perognathus is present on the CP today. Onychomys also occupies arid and semi-arid regions of North America: Onychomys leucogaster, the only species currently found on the CP, prefers shrub steppes and grasslands (McCarty, 1975, 1978). O. leucogaster can be distinguished from its only sister species, Onychomys torridus – found in the Sonoran and part of the Chihuahuan Deserts – using length and width dimension of the m1 and m3 (Carleton and Eshelman, 1979). However, because there is overlap in the measurements for these species (Carlton and Eshelman, 1979), a sample size larger than what is available at ECR2 and RBA is required. Furthermore, many mammal species change their body size in response to climate (Smith and Betancourt, 1998), so morphological criteria are preferable to size for making identifications (Bell et al., 2010).
Taxa with similar relative abundance dynamics group closer together in the NMDS results (Figure 5). Two close taxonomic groups are evident: (1) Dipodomys and Thomomys (highest rank-order abundance across the record, higher relative abundance in the beginning and end of the record, lower in the middle), and (2) Neotoma, Peromyscus, and sciurids (higher in the beginning, then an early dip in abundance followed by an increase, then yet another dip). These groupings do not seem to be connected to climatic affinity because a mixture of taxa with xeric, mesic, or neutral affinity is found in each grouping and abundance does not correlate with MJT.
Desert taxa might be expected to respond more to aridity than to temperature (Stapp, 2010). Several lines of evidence indicate that conditions were relatively cool and wet during time bins A–C and then more arid and warm during D–F (Anderson et al., 2000; Betancourt, 1984; Betancourt and Davis, 1984; Coats et al., 2008; Reheis et al., 2005). Bins E and F encompass a period of time when at least six ‘megadroughts’ took place, including the droughts that are thought to have dispersed the Ancestral Puebloans at Chaco, Mesa Verde, Kayenta, and other large CP settlements (Benson and Berry, 2009; Benson et al., 2007). However, time bins with more similar climate (A–C vs D–F) do not group together on the NMDS, supporting the conclusion that relative and rank-order abundance dynamics are disconnected from climate, at least with regard to the broad-scale trends discussed here. One possible explanation is that time-averaged stratigraphic layers capture short cool or mesic periods within an overall aridifying and warming trend (Hadly, 1996). However, sample sizes for these cool/mesic taxa should be lower than for warm/xeric taxa if taphonomic processes are the same during wet and dry times.
Evenness, as measured using the PIE index, is not significantly correlated with climate. PIE is notably low during bin B (estimated 3569–2700 BP; cool and possibly more wet). At that time, relative abundance of Dipodomys was extremely high, and Neotoma and Peromyscus were at their lowest abundance. Evenness reached a peak in time bin D, when the three most common taxa (Dipodomys, Thomomys, and leporids) all had average relative abundance. Although the uncertainly around PIE does not overlap for most time bins, the magnitude of evenness change is low – from a maximum PIE of 0.85 in D to a minimum of 0.73 in B, a change of only 0.12 – supporting the conclusion that small mammal abundance at these sites was not detectably altered by the changes in climate that occurred during the interval of deposition.
Other factors not tested in this study, like Ancestral Puebloan population size and land use, may influence small mammal diversity dynamics. However, Fisher’s tests reveal that only time bin D was significantly different from the others in terms of the relative abundance distribution across the record. Furthermore, the taxonomic and functional components of this community have not changed markedly through the environmental changes represented by the time sampled by the deposits, suggesting that these groups persist despite natural variation and/or external perturbation in this system. However, finer taxonomic resolution might reveal a signature of climate change: for example, Thomomys subgenus Megascapheus is present in ECR2 and RBA throughout the time period encompassed by this study, but the only specimen of Thomomys subgenus Thomomys is from ECR2 level 3 (time bin F, 644–501 cal. BP). T. subgenus Thomomys is better able to burrow in harder soils than T. subgenus Megascapheus, and soils tend to get harder as climate dries (Marcy et al., 2013). While Thomomys abundance has remained high relative to other taxa for thousands of years, species and subspecies of Thomomys naturally replace one another as climate changes (Hadly and Barnosky, 2009).
Greater sample sizes at RBA and ECR2, particularly at the beginning and end of the record (time bins A and F), could alter the observed abundance patterns and might reveal a correlation between diversity and climate. However, the results reported here are in keeping with other studies of small mammal dynamics through the Holocene in western North America. For example, at Samwell Cave in northern California, Blois et al. (2010) found that small mammal evenness and richness dropped during the Pleistocene–Holocene transition, but during the last ~5 kyr both measures of diversity were fairly stable. At Homestead Cave in northern Utah, Terry and Rowe (2015) showed that energy dynamics, in terms of which small mammal guilds are dominant, changed dramatically at the end of the Pleistocene and again in the last century, but were stable throughout the Holocene. In contrast, Rowe and Terry (2014) compared small mammal diversity dynamics in Homestead Cave and Two Ledges Chamber, a site in northern Nevada, and demonstrated that during the same time interval (8 kyr to the present), species that were shared between sites had entirely different relative abundance dynamics. This suggests that local conditions may have a larger influence on relative abundances at a site than does regional climate (Stegner, 2015).
Conclusion
Two newly excavated late Holocene fossil midden localities, ECR2 and RBA, with 33 AMS 14C dates provide fine-resolution information on small mammal diversity dynamics from ~4.4 to 0.5 kyr in southeastern Utah. These sites comprise a chronology of six time bins. Regional temperature increased across the record, and ‘megadroughts’ marked the time period from ~1100 to 350 BP. I tested for an effect of temperature on community evenness and abundance of 10 small mammal taxa: leporids, perognathines, small sciurids, arvicolines, Cynomys, Neotoma, Dipodomys, Onychomys, Peromyscus, and Thomomys. Despite increases in temperature, measured as reconstructed MJT anomaly, neither is significantly correlated with relative abundance, rank abundance, species richness, or the PIE evenness index. Although abundance does not correlate with climate, certain groups of species have similar diversity trends through time. Identifying what drives these cross-taxa abundance patterns will be important for predicting future abundance changes.
Overall, these results indicate that climatic fluctuations of the magnitude preserved in these deposits do not significantly alter the small mammal community at least on the taxonomic level assessed here. Neither is there evidence that the presence, then exodus, of Native Americans in the region significantly affected small mammals. This lends support to the findings of other studies, like Terry and Rowe (2015), that small mammal communities of western North America were stable throughout the Holocene, but have undergone dramatic changes in structure in the last 100 years as a result of human activities leading to large-scale land cover changes, like invasive annual grasses. The metrics of community structure assessed here can be monitored in the modern small mammal community of this region: if these metrics are significantly different today or in the future from those seen in the fossil record, it may indicate anthropogenic influences outside the range of natural variability.
Footnotes
Acknowledgements
Thanks to AD Barnosky, AR Byrne, EB Davis, EA Hadly, DR Lindberg, and JL Patton for helpful discussion and training. I thank D Erwin, M Goodwin, P Holroyd, and C Marshall of the University of California Museum of Paleontology (UCMP) for curatorial and museum support; T Guilderson and P Zermeno for radiocarbon dating training and assistance; undergraduate assistants M Hardesty-Moore, P Garcialuna, H Bhullar, and P Uppal for processing midden matrix; R Hunt-Foster for assistance with permitting; and the Canyonlands Research Center, Nature Conservancy of Utah, and Redd Family of the Dugout Ranch for support in the field.
Funding
I am grateful to the following sources for funding this research: the Canyonlands Natural History Association, UC Berkeley Department of Integrative Biology, Sigma Xi, Sigma Xi-Berkeley Chapter, the Paleontological Society, the Geological Society of America, the UC Museum of Paleontology, Kappy Wells and the TIDES Foundation, and NSF GRFP Grant DGE-1106400.