Past mass extinction and spatial shifts in biodiversity due to climate change can provide critical data to understand the potential future effects of climate change on global and regional ecosystems. Geologically rapid Quaternary glaciation cycles (Agassiz 1840) are potential analogs for potential biotic change during projected near-future climate change (e.g., Willis and MacDonald 2011; Ordonez 2013). These past changes are believed to have been largely driven by predictable and periodic variation in the earth’s orbit relative to the sun (Milankovitch 1930; Berger 1978; Berger 1980) effecting changes in relative insolation (Imbrie et al. 1992). The global glacial record is well-preserved and continuous over thousands of years in both the Arctic and Antarctic ice-sheets and preserves representative atmospheric gasses that are indicators of past climate (Cuffey & Clow 1997; Alley 2000; Alley 2004). However, North American mid-latitude quantitative climate records are either less continuous or over a short time scale (Viau, Ladd & Gajewski 2012) and primarily limited to ocean and lake sediment analyses including carbonate isotopes and marine diatom assemblages (e.g., Bereiter et al. 2014; Barron et al. 2003; MacDonald et al. 2008), groundwater stable isotope analyses (e.g., Marcott et al. 2013; Lachniet et al. 2014), and scattered biological proxies (e.g., Anderson et al. 2000; Briffa 2000; Cole & Arundel 2005; Wurster et al. 2008). Late Quaternary plant macrofossils from packrat paleomiddens in western North America as presented here, therefore, represent a significant, independent, and novel source of data for large-scale mid-latitude paleoclimate inference for western North America.

North American Packrat Paleomiddens — North American packrat (Neotoma spp.) paleomiddens (fossilized nests) have been a rich source of data for Late Quaternary paleoecological estimates (Betancourt, Van Devender & Martin 1990). Packrats collect both vegetative and reproductive plant material from the local environment to build nests and to establish caches of food. These collections, including leaf and stem fragments, fruits, seeds, and pollen blown into the midden matrix (King & Van Devender 1977) typically form a representative sample (60–90% of taxa) of the local plant diversity within 30–50 m of a midden (Dial & Czaplewski 1990; Finley 1990). Because of the excellent preservation of the materials in dry, protected cavities and under ledges, the fossils can typically be identified morphologically as belonging to extant species from western North America. Paleomidden analysis has long been an important approach to documenting changes in vegetation, population migration and colonization, and morphological and physiological characterization of vegetation of the Late Quaternary as far back as 50,000 years before present (e.g., Wells & Jorgensen 1965; Wells & Berger 1967; Van Devender 1987; Wells & Hunziker 1976; Spaulding & Van Devender 1977; Van Devender 1977; Betancourt & Van Devender 1981; Cole 1982; Cole 1985; Bentancourt et al. 1991; Jackson et al. 2002; Holmgren et al. 2006; Coats et al. 2008; Holmgren et al. 2014; Emslie et al. 2014). Rodent middens preserved in arid environments are often chronologically partitioned, and individual sections can be dated with radioisotope methods to provide a sequential history of vegetation composition (Anderson & Van Devender 1991; Maldonado et al. 2005). The packrat paleomidden fossil record shows remarkably complete biogeographic histories that facilitate the study of past distributional shifts due to local extinction and dispersal processes (e.g., Emslie et al. 2015; Hunter et al. 2001), critical information for understanding how plants respond in natural environments to periods of rapid climate change.

Paleoclimate estimation from Vegetation: Climate Reconstruction Analysis using Coexistence Likelihood Estimation (CRACLE) — The principle that plant species distributions are intimately correlated with climate across broad spatial and temporal scales is universally accepted (e.g., Jackson & Blois 2015). This principle is the basis for a class of climate estimation techniques that use climatic niche overlap between coexisting species to estimate site-specific climate parameters. For example: the Coexistence Approach (Mosbrugger & Utescher 1997) and Mutual Climatic Range (Sinka & Atkinson 1999; Sharpe 2002; Thompson et al. 2012) techniques. In this study, climate parameters are estimated using a maximum likelihood statistical framework for the estimation of climate from vegetation via CRACLE, Climate Reconstruction Analysis using Coexistence Likelihood Estimation (Harbert & Nixon 2015). CRACLE is a relatively new, but a robust method that leverages large online databases of primary biodiversity data such as the Global Biodiversity Information Facility — GBIF ( — and has been shown to provide more accurate estimates of climate than related methods (Harbert & Nixon 2015). Previously, Mutual Climatic Range methods have been applied to the study of Late Quaternary climate for only a few packrat paleomidden samples (Sharpe 2002; Thompson et al. 2012), but a comprehensive regional study has not been completed. This study is the first to apply a taxonomic coexistence climate estimation method to the entire body of data available for Late Quaternary packrat paleomiddens in western North America. It is also the first large-scale application of the CRACLE implementation to paleoclimate modeling (Fletcher et al. 2017).



Estimation of climate was computed based on taxon coexistence and modern distributions via the Climate Reconstruction Analysis Using Coexistence Likelihood Estimation (CRACLE) protocol (Harbert & Nixon 2015). For each taxon in a community, CRACLE calculates a probability function representing the conditional probability of a climate parameter given the presence of that taxon. The joint likelihood function for all co-occurring species is calculated as the sum-log-likelihoods of these individual taxon functions. The maximum of the joint likelihood curve is taken to be the most probable climate value given the association of species and the individual associations with climate. Analyses in this study were done using the ‘vegdistmod’ R package for CRACLE and related methods ( which includes code for CRACLE that is in the process of migrating to the new ‘cRacle’ R package ( The methods used here differ from those originally published (Harbert & Nixon 2015) in two ways. First, a weighted Kernel Density Estimation is used to approximate the conditional probability based on a background data sample. Second, standard bandwidths are calculated using Silverman’s Rule (Silverman 1986).

CRACLE estimates fossil taxon climate tolerance from modern primary biodiversity distribution data using a Smallest Inclusive Group (SIG) approach to identify proxy distributions for extinct or un-identifiable fossil taxa. Here, a SIG is the smallest taxonomic group or clade that a fossil can be reasonably placed in using morphological characteristics. The taxonomic classifications from the USGS North American Packrat Midden Database ( are used as SIGs when the classification is at the genus or species level. One assumption of this approach is that fossils only identified to family level will lack narrow climate signals due to the use of overly broad geographic distributions and likely could be placed into smaller groups. It is expected that taxa in the packrat midden samples represent modern taxa as extinctions during this time are rare, but rearrangements of species into novel communities was common (Emslie et al. 2015; Gill et al. 2009). When generic classifications are used in the CRACLE framework (Harbert & Nixon 2015) the distribution for the entire genus is substituted as a surrogate for the unknown, but extant, species distribution under the assumption of physiological canalization related to climate niche (Huntley et al. 1989; Ackerly 2003; Wake et al. 2009). Using modern species/genus distribution data for fossils that fit within these groups avoids issues of niche shifts and extinction associated with applying this approach to older fossil assemblages that may lack taxa identifiable as equivalent to modern species (Grimm & Potts 2016).

USGS-NOAA Packrat Midden Database—

The work of countless researchers since the discovery of Pleistocene age Neotoma middens (Wells & Jorgensen 1964) to locate, catalog, and identify fossil plant remains in the middens of Neotoma spp. packrats is now largely curated by the United States Geological Survey in the form of the USGS-NOAA North American Packrat Midden Database (available online: The USGS North American Packrat Midden Database ( includes more than 3,000 midden samples, with species identifications, georeferenced localities, and radiocarbon estimated ages of plant material. Because of the historical nature of the dataset, some records lack some information. Sites may have a single age, or more often multiple identified layers from different depositional times that are sampled separately. For our study, the data set was filtered to remove samples without geo-coordinates, with poor coordinate precision, and those with fewer than ten taxa identified to genus or species. The culled data set presented here includes 646 samples (Supplementary material Appendix 2) distributed from the deserts of northern Mexico to the northern Great Basin, eastern California, southern Oregon and Idaho (Supplementary Figure S3). The rarity of the most ancient middens and radiocarbon dating uncertainties restrict usable samples for our study to the last 50,000 years with greater representation of the last 10,000 years (Supplementary Figure S2). On average, there were about 15–20 taxa identified per sample (Supplementary Figure S2), a level of diversity that has been shown to be sufficient for accurate and precise CRACLE models (Harbert & Nixon 2015).

Midden Age Chronology and Paleoclimate Timeline—

Radiocarbon dates for the packrat middens were calibrated using the Northern Hemisphere IntCal13 calibration curve for terrestrial samples (Reimer et al. 2009), and the calibration tools available in the “Bchron” software library (Haslett & Parnell 2008) in R (R Core Team 2016). Age calibration model means and standard deviations were converted to 95% confidence intervals for each site to represent temporal uncertainty in climatic reconstruction. Reconstructions average across all overlapping midden date ranges at each time bin.

Modern climate model and parameters—

Climate data are extracted from the downscaled 2.5-arcminute resolution (~0.041667 degrees) WorldClim (Version 1.4) model grid (Hijmans et al. 2005). WorldClim is a high-resolution continuous grid of interpolated climate data for the world’s land areas derived from a global set of >40,000 weather stations.

An advantage of the CRACLE method for inferring paleoclimate is that the general protocol is flexible and can be applied to a wide variety of climatic parameters. For this study quantitative reconstructions of mean annual temperature (MAT), maximum temperature, minimum temperature, mean precipitation, mean water balance (potential evapotranspiration + mean precipitation), and winter length (number of months with a mean temperature < 5°C) are focused on as these give a broad view of climate change and seasonal variability (Supplementary material Appendix 1).

To better capture nuances of climate-related to drought, the potential evapotranspiration (PET) and water balance (precipitation – PET) were calculated using monthly values for temperature and precipitation from the WorldClim model data (Hijmans et al. 2005) and in R (using code from The Thornthwaite Equation (Thornthwaite 1948; Black 2007) is a suitable model of potential evapotranspiration that relies on data available at the scale of WorldClim including 1) monthly average temperatures in degrees Celcius, 2) mean monthly day-length in hours (calculated from the latitude), and 3) the number of days in each month.

Primary Biodiversity Data

Point-locality distribution data were downloaded from the Global Biodiversity Information Facility API service ( For each taxon (genus or species) queried, localities were removed if any of the following criteria were met: 1) Decimal longitude or latitude was reported as an integer or to only a single decimal place, implying poor precision. 2) The basis of record was listed as “FOSSIL_SPECIMEN”, 3) The locality description had any word matching in-whole or with up to two mismatches “botanic”, “garden”, or “cultivated” to remove observations from cultivated settings. All data were thinned by species to a 20 arcminute grid (the 2.5 arcminute WorldClim grid aggregated by a factor of eight) to reduce the effect of spatially biased collection patterns. Taxa with fewer than five records after filtering were excluded as incomplete data. The final occurrence dataset consisted of 368,702 unique records for the 341 genera and 533 species observed in the 646 studied midden samples.

Regional CRACLE Validation Testing

To test whether CRACLE is adequate for capturing climate information at a regional scale comparable to the total climate space for the Late Quaternary of western North America the model was applied to vegetation surveys from the United States, Canada, and Mexico mobilized by the Botanical Information and Ecology Network (BIEN, via the R library ‘rbien’ (Maitner et al. 2018). For each survey, primary biodiversity data from GBIF was downloaded as described above and all the same filtering and analytical steps were applied as with the USGS Packrat Midden Database. In the final distribution data, there were 5,563 species and 915,414 records. One additional filter was applied for the modern sites which was to remove surveys with reported elevation greater than 200 meters above or below the estimated elevation used in the WorldClim model (Hijmans et al. 2005; Elevation filtering will reduce calculated anomalies by removing surveys that are either poorly georeferenced or were collected on a very strong elevation gradient that is not captured well in the WorldClim climate interpolation. All code and data to reproduce this validation test are publicly available (

Results and Discussion

CRACLE Modern Validation Experiment — To validate CRACLE as a sensitive and robust method for estimating climate across the Late Quaternary of western North America a modern set of more than 3,493 vegetation plot surveys covering North America west of 100°W and between 25°N and 50°N (mobilized via the Botanical Information and Ecology Network — BIEN — database were analyzed using the same CRACLE methodology as for the packrat midden fossil communities.

These modern sites cover between –3°C to 23°C degrees of mean annual temperature and 100–1,100 mm of annual rainfall, a latitudinal gradient of ~20 degrees, and altitudes between 170 and 4000 meters. This dataset should well represent the scope of climate space for the Late Quaternary of this part of western North America.

Temperature variables (e.g., Mean Annual Temperature — MAT) were estimated with high correlation (ρ = 0.93) and low error rates (median absolute anomaly = 1.3°C — Figure S1; Supplementary Table S1) relative to the WorldClim estimates for those sites. Estimation of precipitation parameters (mean annual precipitation and annual water balance) exhibited greater variation in results, but with median anomalies of 79 mm for mean annual precipitation and 96 mm for annual water balance (Table S1). Seasonal estimation, represented in this study by winter length (the number of months with mean temperatures below 5°C), performs well with a median anomaly of only 0.05 months, or about 1.5 days. The performance in this regional dataset is comparable to the global set of test sites originally published for CRACLE (Harbert & Nixon 2015) but from a dataset almost 25 times larger. These results show that CRACLE is a robust tool in for this regional flora and the climate space relevant to the packrat midden macrofossil record of the Late Quaternary.

Paleoclimate Overview — The climate record inferred by CRACLE analysis of the 646 dated western North American packrat paleomidden samples is consistent with patterns of northern hemisphere temperature changes over the last ~50,000 years (Figure 1A, B) based on other methods and proxy data (Cuffey & Clow 1997; Alley 2000; Alley 2004; Shakun et al. 2012; Marcott et al. 2013). Coldest temperatures with anomalies of –2 to –4°C were estimated between 22 ka and 15 ka, during and just after the last glacial maximum (LGM – 26–19 ka; Clark et al. 2009) in Western North America (Figure 1A).

Figure 1 

50,000 years of CRACLE Mean Annual Temperature record, Northern Hemisphere Ice-Core estimated temperature data, peak summer insolation at 65°N, and atmospheric CO2 concentration. A) Sliding window 1000 year moving average Mean Annual Temperature (MAT) difference estimated by CRACLE is shown as the mean deviation +/–2 standard error (dashed lines) from modern average for the 646 study sites. B) Northern Hemisphere multiproxy temperature anomaly stack for the Late Pleistocene (Shakun et al. 2012) and Holocene (Marcott et al. 2013). C) Greenland air temperature sliding window 1000 year moving average estimated anomaly (Data from the Greenland Ice Sheet Project 2, & Clow 1997; Alley 2000; Alley 2004). D) Mean daily summer (June) insolation at 65°N (kJ/m2) estimated following the geometric method (Willis & MacDonald 2011) using the ‘palinsol’ R library (Crucifix 2016). E) Mean CO2 concentration (ppm) estimated from a composite of Antarctic ice core records (Alley 2000).

The late Pleistocene temperature minima were followed by a general warming trend that continued until after the onset of the Holocene likely driven by the collapse of the continental ice-sheets and peak northern latitude summer insolation (Figure 1D). Peak positive anomalies of +1.5°C above modern values consistent with the Holocene Thermal Optimum (Renssen et al. 2012) at 7–8 ka (Figure 1A).

Atmospheric CO2 recorded in the Antarctic ice cores (Figure 1E) (Bereiter et al. 2015) increases ahead of mean annual temperature increases in the CRACLE inferred record (Figure 1A) supporting previous observations of global climate change at the end of the Pleistocene (Shakun et al. 2012). Although showing dramatic extremes (~6°C total increase in temperature), the magnitude of CRACLE-packrat estimated mean temperature shifts at these latitudes is much less than is observed in the reliable ice-core proxies of Greenland (Figure 1A, C) (Cuffey & Clow 1997; Alley 2000; Alley 2004). A pattern that is consistent with a scientific consensus that global temperature change and ecosystem response is typically greater in polar regions during past warming events (Renssen et al. 2012; Kiessling et al. 2012).

The rate of Climate Change —The maximum rate of mean annual temperature change estimated in this study (~6°C/6 to 8 thousand years = 0.075–0.1°C/100 years; Figure 1) is a slower rate than current projections of near-future warming scenarios for western North America. The projected rate of warming across the western United States through the year 2100 is projected to be between 1 and 4°C in all IPCC Representative Concentration Pathways (Pachauri et al. 2014 – see page 61), including the stringent mitigation scenario —RCP 2.6— that is, at this point, very unlikely to be attained (Pachauri et al. 2014). In the Southwestern United States, recent historical records suggest increasing temperatures are now occurring at a rate of ~3–4°C/100 years, a rate consistent with or exceeding model projections for the next century (Cole 2016). Based on these conservative estimates contrasted against CRACLE paleomidden estimates, projected warming in the next century is likely to occur at a rate 10–20? greater than any millennial average in the past 50,000 years (Cole 2016). This disconnect highlights the need for studies that infer potential future biotic change (e.g., Mathys et al. 2017) and migration rates (Roberts & Hamman 2016) to consider the unparalleled effects of very rapid periods of climatic change rather than just the total magnitude of projected change (Cole 2010).

Pleistocene Pluvial Lakes — One of the most striking geologic features of the late Pleistocene is the occurrence of large lake systems in the Great Basin during the Last Glacial Maximum (Reheis 1999; Lyle et al. 2012; Munroe & Laabs 2013; Reheis et al. 2014). Estimates of mean annual precipitation (MAP) and water balance (MAP – potential evaporation) for the LGM are consistent with a wet-glacial hypothesis (Oster et al. 2015), with additional rainfall of 50–100 mm/year over modern values (Figure 2C). Driven by cooler temperatures and reduced evaporation, relative water balance (i.e., subtracting estimated evaporation) was greater than modern values by +200 mm/year to +400 mm/year from 22 ka to at least 15 ka (Figure 2D). During the Bolling-Allerod interstadial (a pronounced Northern Hemisphere warming episode from 14.7 ka to 12.7 ka), the wetter regime across the region started to collapse. At this time, CRACLE paleomidden analyses show rapidly diminishing rainfall from 15 ka towards the modern regime by ~12 ka. This decrease in rainfall is followed closely, though with less abrupt changes, by water balance decreases as temperatures increased towards the HTO (Figure 2A, 2B). Interestingly, CRACLE-estimated winter length remained longer than the modern by about one month until around 9 ka (Figure 2E), despite rising annual (Figure 1A) and seasonal extreme temperatures (Figure 2A, 2B).

Figure 2 

Temperature extremes, precipitation, water balance, and winter season CRACLE anomaly time series. A) Maximum annual temperature (°C), B) minimum annual temperature (°C), C) mean annual precipitation (mm), D) mean annual water balance (mm) – potential evapotranspiration + precipitation, E) estimated winter length (number of months with mean temperature less than 5°C). Solid lines indicate mean difference; dashed lines are +/–2*SE. Shaded regions indicate time periods where significant increases (blue) or decreases (purple) in taxonomic diversity are observed (Supplementary Figure S1).

The rise of the Pleistocene pluvial lakes in the Great Basin has been attributed to an influx of moisture due to increased rainfall and/or limited evaporation by the full glacial cold and reduced summer insolation (MacDonald et al. 2008; Lyle et al. 2012; Ibarra et al. 2014; Lachniet et al. 2014; Oster et al. 2015). Our estimated changes in water balance during the LGM (22 to 15 ka) show a positive linear relationship with latitude — larger positive anomalies in water balance to the north — (Figure 3). The positive trend with latitude is consistent with other contemporary models that show greater rates of temperature increase in northern, more polar latitudes than in southern or more tropical latitudes (Renssen et al. 2012; Kiessling et al. 2012). Predicted differences in water balance were used to adjust the modern climate model grid using a linear model of latitude vs. estimated water balance change (Figure 3) for the samples from the LGM maximum positive anomalies and mid-Holocene maximum negative anomalies. During the LGM these models produce anomalies of up to +500 mm/year of water balance, with greater values in the north, pushing positive water balance regimes to the south and down-slope.

Figure 3 

Linear modeling of Great Basin water balance (PET + MAP) during the estimated periods of greatest temperature departure from modern. A positive relationship exists between CRACLE inferred water balance anomalies vs. latitude for both the temperature minimum of the from 22 ka to 15 ka during the Last Glacial Maximum (r2 = 0.0927; p = 0.004109; n = 77) and the Holocene thermal optima between 8 ky and 2 ka (r2 = 0.1053; p = 0.000; n = 176) exhibit a positive and near-linear relationship. Applying the linear models to predict the water balance across the study region yields continuous adjusted estimates of water balance for the LGM and Holocene and visual comparison to the modern (1950–2000) average (Hijmans et al. 2005). Maps were generated with R version 3.5.1 (R Core Team 2018) and associated geospatial libraries raster version 2.6-7 (Hijmans 2017) and rasterVis version 0.41 (Lamigueiro & Hijmans 2018).

The Holocene: Implications for future drought scenarios — Based on CRACLE estimates, much of the Holocene is characterized as warmer by 1–2°C, and drier than the present (Figure 1A). Our reconstructions of precipitation and water balance suggest that the Holocene after ~8 ka was as dry as, or dryer than, the current climate throughout the western U.S. However, it is important to note that current near-future climate projections indicate we will exceed this worst-case historical scenario, and perhaps more importantly, at a much quicker rate. Temperature maxima in the Holocene are estimated to occur at 7–8 ka and 2–4 ka for all temperature parameters (Figure 2A, B). In turn, temperature increases appear to drive estimated water balance minima at those same time periods (Figure 2D). This result is in contrast to initial warming followed by a steady downward trend in MAT between the HTO and the present reported from multiproxy climate reconstruction syntheses (Figure 1B; Marcott et al. 2013). The CRACLE estimated warm periods and the pattern of associated expansion of drought regimes (Figure 3) are reasonable analogs for future warm climates in arid and semi-arid regions of western North America. CRACLE estimates of drought (negative water balance) suggest that parts of the southern half of our study area, including Southern California, Arizona, and northern Mexico, may see shifts from positive to negative water balance regimes if the projected warming by the year 2100 verifies (Figure 3).

Estimation of Holocene water balance also reveals that these warm and dry climate periods had less impact to the north of the study area. In particular, the northern Sierra Nevada range and coastal regions of northern California, Oregon, and Washington show little change in the overall wet climate through the Holocene (Figure 3). This moisture regime relative stability may imply that these regions are historical, and potential future refugia. It is interesting that these regions (such as the Siskiyou Range) are considered refugia for many plants that were more widespread in the past (Steele & Storfer 2006; Ackerly et al. 2010; Olson et al. 2012).

Holocene Temperature Variability — All estimates of temperature presented here (mean annual, minimum, and maximum) show marked variability during the Holocene (Figures 1A, 2A, B) including a local temperature minimum around 4–5 ka. The dominant global temperature trend of the Holocene has been temperature stability except for the last 2000 years where anomalous cooling has been inferred (Marcott et al. 2013; Marsicek et al. 2018). However, climate variability at local and regional scales can get lost in such broad syntheses. CRACLE analysis reconstructs an overall warm Holocene with a cooling trend over the last 2,000 years of the Holocene consistent with these global trends (Marcott et al. 2013).

The CRACLE temperature estimates show first a ~2°C cooling trend between 5–7 ka followed by an equal warming trend between 3–5 ka (Figures 1A, 2A, B). Model data suggest a 2–4°C maximum temperature anomaly during the early Holocene (Renssen et al. 2012) in this region, an inference corroborated by the early warm period in the CRACLE estimates (Figure 1A). Sea surface temperature (SST) along the Northern California coast may have been cooler than present during the mid-Holocene (Barron et al. 2003). Diatom records from the northern Gulf of Mexico suggest an amplified North American Monsoon (NAM) resulting in a wet summer pattern during the mid-Holocene as well (Poore et al. 2005). Cooler SST (Barron et al. 2003) and increased summer rainfall (Poore et al. 2005) may help explain the temperature minimum at 5 ka in our results as a regional climatic response.

Millennial Climate Variation and Plant Biodiversity Dynamics — Our CRACLE estimates for the end of the Pleistocene show that mean annual temperature increased from the LGM minimum prior to decreases in precipitation near the Pleistocene-Holocene transition (Figure 1). Surprisingly, this period also displays a significant millennial-scale increase in average within-midden plant diversity —temporal α-diversity— (data and statistical test described in Supplementary Figure S2). Plant diversity in the arid regions of the Southwest is lowest in dense, high elevation forest, low elevation grasslands, and desert scrub, but highest in open, dry woodlands due to a mixing of both woodland and grassland/desert elements in this zone (Whittaker & Neiring 1975). At the end of the Pleistocene, a warming trend with moderate rainfall that persisted across large portions of this region could have driven an expansion of a semiarid climate with open woodland vegetation types. Thus, increases in open woodland vegetation with relatively higher species diversity is one potential mechanism for the observed increase in within-midden diversity in this period (Supplementary Figure S2).


The Late Quaternary plant macrofossil record preserved in packrat paleomiddens provides exceptionally strong evidence to estimate a detailed and well-supported timeline of climatic change across western North America. This study is the first application of the CRACLE approach to Quaternary paleoclimate estimation and the first large-scale quantitative reconstruction of climate in this region using plant macrofossils and demonstrates the potential application of CRACLE to other paleoecological records and databases (e.g., the Neotoma Paleoecology Database — Williams et al. 2018). The 50,000-year climate profile estimated using CRACLE includes many well-understood features of global climate history found with other methods and in other regions, validating the overall approach and its application in this region. Perhaps more importantly, this reconstruction quantifies changes in numerous climate variables not typically accessible through other proxy methods (e.g., water balance and winter length). The climate profile also sheds light on key events of the geologic and fossil record including the development of large pluvial lakes in the Great Basin during the terminal Pleistocene glacial cycle (>17 ka) during which we reconstruct increases in overall precipitation and reduced evaporation. The rapid temperature increase, a decrease in precipitation, reduced water balance and implied greater drought marking the Pleistocene-Holocene transition (12 ka–9 ka) is evidence of a period of rapid climate change that elsewhere has been implicated in coordinated extinctions of Pleistocene megafauna in North America (Guthrie 2006; Zazula et al. 2014). However, future projected anthropogenic climate change might occur at a rate of ca. 1–4°C per century (Pachauri et al. 2014; Cole 2016), at least one order of magnitude greater than the millennial rate of about 0.1°C per century estimated here for the Pleistocene-Holocene transition. The rate estimated in this study for the Pleistocene-Holocene transition is less than the ca. 0.5°C per century for polar regions for the same period inferred from the GISP2 cores (Cuffey & Clow 1997; Alley 2000; Alley 2004). This rate disparity highlights the utility of the continuous mid-latitude climate reconstructions made by CRACLE using the packrat midden macrofossil system (or other similar biological data sources) for aiding the quantification of the latitudinal climate change differential. Climate shifts in the past have resulted mostly in rearrangements of the standing biodiversity of the regional flora of western North America rather than extinction. The generally modest rates of climatic change may have allowed for processes of migration and dispersal to track environmental shifts over long periods. Given the rapid changes projected for the next century, there appears to be no broad analog in the Late Quaternary of western North America for the for the projected >1°C per century temperature increase that the world faces in the next century.

Additional Files

The additional files for this article can be found as follows:

Figure S1

Modern validation CRACLE Anomalies. DOI:

Figure S2

Midden count and mean taxon diversity timeline. DOI:

Figure S3

Geographic distribution of 474 packrat midden fossil localities used for CRACLE paleoclimate proxy. DOI:

Table S1

Modern CRACLE performance statistics. DOI:

Appendix 1

Results file for MAT, MinT, MaxT, MAP, wbalann, winter length by site. DOI:

Appendix 2

Sites used (Localities + links to USGS database). DOI: