Spatial sampling bias in the Neotoma paleoecological archives affects species paleo-distribution models

Quaternary Science Reviews 198 (2018) 115e125

Contents lists available at ScienceDirect

Quaternary Science Reviews journal homepage: www.elsevier.com/locate/quascirev

Spatial sampling bias in the Neotoma paleoecological archives affects species paleo-distribution models Richard Inman a, b, *, Janet Franklin c, Todd Esque a, Kenneth Nussear d a

US Geological Survey, Western Ecological Research Center, Las Vegas Field Station, 160 N. Stephanie Street, Henderson, NV, 89074, USA School of Geographical Sciences and Urban Planning, Arizona State University, PO Box 875302, Tempe, AZ, 85287, USA c Department of Botany and Plant Sciences, University of California e Riverside, 900 University Ave, Riverside, CA, 92521, USA d Department of Geography, University of Nevada e Reno, 1664 N. Virginia Street, Reno, NV, 89557, USA b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 14 May 2018 Received in revised form 15 August 2018 Accepted 16 August 2018

The ability to infer paleo-distributions with limited knowledge of absence makes species distribution modeling (SDM) a useful tool for exploring paleobiogeographic questions. Spatial sampling bias is a known issue when modeling extant species. Here we quantify the spatial sampling bias in a North American packrat midden archive and explore its impact on estimating paleo-distributions. We test whether (1) spatial sampling bias inherent in this macrofossil record can influence estimates of paleodistributions, (2) this bias can alter the ability to measure shifts in distributions and climatic niche breadth from the Northgrippian subdivision of the Holocene (8.3 ka e 4.2 ka) to present day (1950 e2000 yr), and (3) bias correction methods can improve estimates of paleo-distributions and analyses of range shifts and niche breadth. We estimate spatial sampling bias for the mid-Holocene period with a three-stage statistical model, each representing a hypothesized source of bias: fossil site availability, preservation and accessibility. This approach enables the use of SDM to evaluate three separate paleodistributions calibrated on the packrat midden archive: those without bias correction (s-naïve), those created with a standard method (s-standard), and those created with a novel alternative (s-modeled) incorporating the three-stage model of bias. We find that paleo-distributions modeled for the midHolocene without bias correction (s-naïve) provided poor estimates of hindcast paleo-distributions, and that the s-modeled correction method improved paleo-distributions for our six species with, on average, 50% higher overlap to hindcast distributions than s-naïve paleo-distributions (s-standard results fell between s-naïve and s-modeled). Published by Elsevier Ltd.

Keywords: Holocene Paleogeography North America Data treatment Data analysis Neotoma Spatial sampling bias Species distribution modeling Niche breadth Macrofossils

1. Introduction A core focus of biogeography rests in understanding the determinants of species distributions and the processes by which they change. Towards that goal, rapid development of paleoecological archives and analytical tools over the past 20 years has enabled investigations of broad macro-ecological, evolutionary and conservation questions about the mechanisms and forces altering patterns of biodiversity throughout the history of our planet (Brewer et al., 2012; Swetnam and Allen, 1999). Species distribution modeling (SDM) has become a widely used tool in paleobiogeographic studies because in a presence-background

* Corresponding author. US Geological Survey, Western Ecological Research Center, Las Vegas Field Station, 160 N. Stephanie Street, Henderson, NV, 89074, USA. E-mail address: [email protected] (R. Inman). https://doi.org/10.1016/j.quascirev.2018.08.015 0277-3791/Published by Elsevier Ltd.

modeling framework, SDM does not rely on explicit knowledge of absence localities (Franklin, 2010), which are often difficult to determine (Brotons et al., 2004). The ability to use presence-only data, such as fossil records, allows SDM to infer paleodistributions with limited knowledge of absence and address paleobiogeographic questions of niche stability (e.g. Stigall, 2012), s-Bravo et al., 2008; Veloz et al., 2012), range dynamics (e.g. Nogue speciation (e.g. Peterson and Ny ari, 2008), and extinction (e.g. Lorenzen et al., 2011), among many others. When used to investigate paleo-distributions, SDM offers two analytical opportunities: 1) validation of models calibrated on extant species that have been ‘hindcast’ to past environmental conditions, and 2) direct calibration with past environmental conditions (paleo-SDM). The former has been used most often, driven by the wealth of spatially explicit ecological archives available for extant species (Moreno-Amat et al., 2017). In this approach, models

116

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125

are calibrated using present-day conditions and are applied to environmental data representing paleoclimatic conditions to create spatial predictions for the historic period of interest. Paleoecological archives are then used to quantify the accuracy of these hinds-Bravo, 2009) or cast projections (e.g. Franklin et al., 2015; Nogue to qualitatively evaluate the congruence between hindcast projections and fossil localities (e.g. Carnaval and Moritz, 2008). We refer to these predictions of paleo-distributions as ‘hindcast’ distributions. These hindcast projections describe areas where a species could have occurred given its present-day realized niche; but rarely do these hindcast distributions consider exogenous factors such as fire or megafaunal disturbances that may have constrained s-Bravo, 2009). In paleo-distributions (Gill et al., 2009; Nogue contrast, paleo-SDM uses paleoecological archives for calibrating models under past environmental conditions. Paleo-SDM is used less often, likely because paleo-distribution data are usually sparse in time and space and can be poorly resolved chronologically and taxonomically (Moreno-Amat et al., 2017). However, paleo-SDM allows a key assumption of hindcasting to be relaxed, namely, the requisite of niche conservation through time. Under niche conservation, species-environment parameters are maintained through time, even with environmental change (Wiens and Graham, 2005). By directly calibrating on fossil archives, paleo-SDM can estimate new species-environment parameters for different chronological periods. Each unique set of parameters may provide evidence of a changed niche, and in conjunction with spatial predictions of habitat potential, can be compared across different time slices to evaluate potential change in niche and habitat (e.g. Saupe et al., 2014). Often, studies using paleo-SDM rely on the spatial predictions of habitat potential because many clear metrics exist for evaluating overlap among time slices, such as the Sorenson's similarity index (Sørensen, 1948), Schoener's D (Schoener, 1968) or Godsoe's ESP (Godsoe, 2013). These metrics, among others, allow for simple comparisons of spatial predictions rather than complex assessments of niche complexity, dimensionality, or breadth. When coupled with hindcasting, direct tests of niche conservation can be made between distributions derived from hindcasting and those from fossil calibrated models (e.g. Veloz et al., 2012; Walls and Stigall, 2011). These coupled approaches have been used often, contributing evidence that niches evolve slowly (e.g. MartinezMeyer et al., 2004; Peterson, 2011; Stigall, 2012). However, paleoecological archives used for paleo-SDM often have small sample sizes, which can be biased in space and time (Varela et al., 2011; Vilhena and Smith, 2013). Spatial bias is problematic for paleo-SDM because these methods assume niches are sampled over the full range of environmental conditions in which they occur (Phillips et al., 2009) and that sampling mirrors the background distribution of environmental covariates (Araújo and Guisan, 2006). These assumptions are not often met with paleoecological archives due to spatial variation of taphonomic conditions in different deposits (Allison and Bottjer, 2010), which results in spatially biased fossils where more fossils are found in certain areas due not to a greater prevalence of an organism, but instead due to a lack of fossils in other areas. Spatial bias in biodiversity data for extant species is often addressed with one or more bias correction methods that have been developed for presence-background frameworks. These methods stem from careful filtering of observation data (e.g. Boria et al., 2014) or from estimating the biased sampling distribution (s) and manipulating background selection weights to result in proportional background samples (e.g. Phillips et al., 2009). In order be effective, these bias correction methods require large sample sizes. Paleoecological archives, however, rarely offer large samples (Moreno-Amat et al., 2017; Varela et al., 2011), and the ability to estimate s from them is not often tested. An alternative method to

estimate s that does not rely on large sample sizes may improve estimates of paleo-distributions, thereby allowing analyses of range shifts and niche characteristics through time. Here we explore potential effects of sampling bias on paleo-SDM and investigate an alternative for estimating s with a focal group of extant plant species. Specifically, we test whether (1) spatial sampling bias inherent in a commonly used paleoecological archive, the USGS/NOAA North American Packrat Midden database (here after NAPM database; Strickland et al., 2013), can influence estimates of paleo-distributions, (2) this bias can alter our ability to measure shifts in distributions and niche breadth from the Northgrippian subdivision of the Holocene (8.3 ka e 4.2 ka) (here after midHolocene; Cohen et al., 2013) to present day (1950e2000 yr), and (3) bias correction methods can improve paleo-distributions and analyses of range shifts and niche breadth. While there is evidence that the niches of many North American plant species have shifted since the last glacial period of the late Quaternary, changes since the mid-Holocene period have been minimal (Ordonez, 2013; Veloz et al., 2012). We therefore use paleo-distributions created by hindcasting present-day models to the mid-Holocene as our reference distributions and compare them to three separate paleodistributions calibrated on the fossil record: those without bias correction (s-naïve), those created with a standard method (sstandard), and those created with a novel alternative (s-modeled). Using hindcast distributions for paleo-SDM evaluation assumes that hindcast distributions are representative of mid-Holocene distributions, and so to independently verify the paleo-SDM distributions we also used pollen records from lake sediment cores from the western USA not used to calibrate paleo-distributions. We aim to identify an effective method for reducing sampling bias in paleo-SDM and to highlight how bias may affect analyses of range shifts and niche breadth. 2. Materials and methods 2.1. Study area and environmental gradients Our study region covered 3,171,335 km2 of the western USA encompassing the locations of packrat middens represented in the NAPM database (Strickland et al., 2013). We assembled uncorrelated raster data describing 11 climatic and physiographic environmental conditions to characterize present-day (1950e2000 yr) and mid-Holocene time periods, and generalized each at a spatial scale of 1 km (Appendix A). Interpolated climate variables for current conditions were obtained from WorldClim (Hijmans et al., 2005) as average monthly values and climate variables representing the mid-Holocene were obtained as downscaled paleoclimate simulations from the National Center for Atmospheric Research (NCAR) from the Community Climate System Model (CCSM4; Gent et al., 2011), also available through WorldClim. We assumed that most changes between the two periods were limited to climatic variables and that any changes in surface physiographic conditions due to Holocene erosion or deposition were minimal and within 2x the maximum vertical error of the terrain elevation data; 10 m (Danielson and Gesch, 2011). 2.2. Species distribution modeling We developed present-day and paleo-distributions for six extant species well represented in the NAPM database that occur in a range of habitat types, including two small perennial shrubs inhabiting cold desert habitats (Artemisia tridentata, Coleogyne ramosissima), a small deciduous tree inhabiting foothills and low mountain elevations (Quercus gambelli), a small conifer tree confined to northern latitudes (Juniperus communis), and two large

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125

conifers occurring in high elevation mountain habitat (Abies concolor, Pinus ponderosa). The NAPM database has been used extensively to assess regional vegetation community shifts (e.g. Jackson et al., 2005; Thompson and Anderson, 2000), paleoclimate (e.g. Arundel, 2002; Jacobson, 1991), and faunal community composition (e.g. Mead, 1981; Van Devender and Mead, 1978) during the late Pleistocene and throughout the Holocene. More recently, packrat midden archives have been used to estimate spatially explicit paleo-distributions for some species (e.g. Angulo et al., 2017), though few studies have used paleo-SDM methods. The paucity of studies using paleo-SDM is likely due to the spatial bias inherent in the packrat midden archive (Mensing et al., 2000; Webb and Betancourt, 1990). We develop a novel approach to estimate and correct for this bias (s) so that paleo-SDM may be used to create spatially explicit paleo-distributions from the NAPM database. We obtained macrofossil observations from the NAPM database (v.4; Strickland et al., 2013) and calibrated 14C ages into calendar ages using the Bchron package (Haslett and Parnell, 2008) in R 3.3.3 (R Core Team, 2016) and selected records from the mid-Holocene period to create paleo-distributions with each of the three paleoSDM methods. A species was considered present at a given midden location if any of the calibrated ages associated with macrofossils for that species were dated to the mid-Holocene, resulting in sample sizes ranging from 4 to 22 (Appendix B). Extant locations from 1950 to 2000 for the six species were obtained from the Geographic Biodiversity Information Facility (GBIF, www.GBIF.org; Appendix C) in October 2017, resulting in sample sizes ranging from 376 to 2039 (Appendix B). 2.3. Present-day distributions We used MaxEnt version 3.4.0 (Phillips et al., 2006), a widelyused software for SDM in presence-background frameworks, to create distribution models for present-day conditions. For each species, we created a model with all 11 environmental explanatory variables and sequentially removed those contributing the least to model fit using a step-wise jackknife test of training gain (Elith et al., 2011). We stopped removing variables when a noticeable drop in the Area Under the receiver operating characteristic Curve (AUC; Fielding and Bell, 1997) was observed with 20% withheld test data. The resulting present-day model for each species contained three to five explanatory variables. We used the complementary log-log function and allowed inclusion of each feature class (linear, quadratic, product, threshold and hinge) in a bootstrap framework with 100 iterations, and saved the standard deviation across all iterations to approximate model error at each grid cell. We use the FactorBiasOut (hereafter, s-standard) algorithm (Phillips et al., 2009) for reducing spatial sampling bias inherent in present-day GBIF observations. The s-standard method estimates the species distribution (p) by deriving the combined distribution of p and biased sampling distribution (s), and then factoring s out. This method relies on knowledge of s, which, in practice, is rarely known. However, because the biased observation dataset is a sample of s (i.e. observations are sampled from the biased sampling distribution), s can be estimated when the observation dataset is large (Phillips et al., 2009). In these cases, s is represented as an auxiliary variable, often in the form of a sampling intensity bias grid (Elith et al., 2011). We develop a bias grid to estimate s by creating a kernel density raster from the biased observation dataset of each species, and estimate the bandwidth for each kernel with a crossvalidated selection method to minimize mean-square error (Baddeley et al., 2015). The resulting bias grid was rescaled to 1 to 20, to give greater selection probability to areas with higher densities of observations as recommended by Elith et al. (2011).

117

2.4. Mid-holocene distributions Paleo-distributions were created with 4 methods: 1) hindcasting from present-day conditions, 2) paleo-SDM without bias correction (s-naïve), 3) paleo-SDM with standard correction (sstandard), and 4) paleo-SDM with model correction (s-modeled). Hindcast paleo-distributions were created by projecting presentday distributions with the same explanatory variables, but with climate values from the NCAR CCSM4 paleoclimate simulations representing the mid-Holocene. In order to validate each hindcast paleo-distribution, we calculate a test AUC score using the macrofossils obtained from the NAPM database since they were not used in model calibration. In contrast, the three paleo-SDM methods relied on calibrating new models with the macrofossil records obtained from the NAPM database dated to the mid-Holocene period. For each of the paleo-SDM methods, we used the same bootstrap framework in MaxEnt with 100 iterations and saved the standard deviation across all iterations to approximate model error. Due to the small sample sizes of macrofossil records, we did not withhold 20% test data as with present-day models. The s-naïve paleo-distributions were created without any bias correction, while the s-standard used the FactorBiasOut algorithm (described above) with a kernel density bias grid estimated from the macrofossils to approximate the biased sampling distribution of the NAPM database. Because macrofossil sample sizes can be small, we expected that the s-standard method would not sufficiently reduce spatial sampling bias and therefore developed an alternative approach, smodeled, to reduce model bias. We hypothesized that the macrofossils obtained from the NAPM database would be biased towards areas that 1) were suitable for Neotoma (packrat) populations during the mid-Holocene, and 2) have been suitable for fossil preservation since the mid-Holocene. We therefore develop a new approach to reduce spatial sampling bias by estimating the biased sampling distribution (s), and using it to manipulate background samples in MaxEnt. 2.5. A new approach to paleo-SDM bias correction: s-modeled Our s-modeled approach draws on three separate statistical models to estimate s, each representing a hypothesized source of bias in fossil data: availability, preservation and accessibility. The first model, availability, accounts for the prerequisite of Neotoma species being present at a midden location during the midHolocene. This component was modeled using extant observations obtained from GBIF for six species of Neotoma (N. albigula, N. cirerea, N. lepida, N. leucodon, N. mexicana and N. micropus) to create a present-day distribution for the Neotoma species assemblage with MaxEnt version 3.4.0 (Phillips et al., 2006). As with other implementations of Maxent in this study, we use a bootstrap framework with 100 iterations and 20% of sample observations withheld for testing with AUC. Environmental explanatory variables consisted of the same datasets described in Appendix A, but were extended to include the conterminous USA to include the full range of Neotoma. Variable selection was conducted with a jackknife test of gain. Biased sampling distributions were accounted for using the s-standard bias correction method because the sampling distribution is suitable to support this approach; sample size of extant observations used to calibrate present-day models of Neotoma was 6642. A paleo-distribution for the mid-Holocene was created by projecting the present-day distribution onto climate values representing the mid-Holocene (hindcasting). We validated these hindcast paleo-distributions using AUC with 202 macrofossil locations of Neotoma obtained from FAUNMAP (FAUNMAP Working Group, 1994) dated to the mid-Holocene, and provide a summary of

118

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125

variable contributions in Appendix D. The second component, preservation, accounts for the physical conditions needed to create a well-preserved midden, which are limited to physiographic and substrate conditions such as rocky hillsides, cliff faces, or talus slopes, and are usually on north, northeast, south or southwest facing slopes (Webb and Betancourt, 1990). Arid conditions are also needed for midden preservation, such that locations with high soil moisture or dense vegetation canopy cover are not very suitable (Webb and Betancourt, 1990). We considered explanatory variables representing climate (e.g. aridity and temperature), physiography (e.g. solar exposure, drainage, ruggedness, slope and aspect), underlying material composition (e.g. consolidated and unconsolidated), and geologic characteristics (e.g. sedimentary-, metamorphic-, and igneous-rock types) for the conterminous USA (Appendix E). We calibrated this model on locations of middens containing records of any plant species dated to the mid-Holocene. Again, we use a bootstrap framework with 100 iterations and 20% of observations withheld for testing with AUC, and selected variables using a jackknife test of gain. We also provide a summary of variable contributions in Appendix D. The third component, discovery, represented the potential for macrofossils to be discovered, and was therefore calibrated on all known midden sites in North America. We considered explanatory variables hypothesized to influence fossil discovery, such as erosional proxies, human accessibility and vegetation surface conditions (Appendix F). Here again, we use a bootstrap framework with 100 iterations and 20% of observations withheld for testing with AUC and select variables using a jackknife test of gain. A summary of variable contributions is provided in Appendix D. The product of the three statistical models for availability, preservation and discovery components was used to represent s in the smodeled paleo-distributions for each species and was used to select background records with a probability equal to s.

2.6. Paleo-SDM validation with hindcast distributions: permutation tests We hypothesized that the s-modeled paleo-distributions would show higher congruence to the hindcast paleo-distributions than would either the s-naïve or s-standard methods. This hypothesis relies on the assumption that hindcast distributions, calibrated on large samples from present-day conditions, should represent true paleo-distributions better than the paleo-SDM models (which are calibrated on biased and small sample sizes). It is reasonable to assume that very little niche evolution occurred since the midHolocene (Peterson, 2011; Stigall, 2012), especially when a species' niche is characterized broadly in climatic and environmental space (Wiens et al., 2009). Therefore, hindcast distributions over short intervals present opportunities for assessing the effectiveness of bias correction methods. To assess congruences between hindcast and the s-naïve, sstandard and s-modeled paleo-distributions, we used the Expected fraction of Shared Presences overlap metric (ESP; Godsoe, 2013) by comparing the estimated habitat potential of each across all cells with the equation:

P 2 j P1j P2j  ESP ¼ P  j P1j P2j

(1)

where P1j denotes the habitat potential for the hindcast paleodistribution at location j, and P2j denotes the prediction generated from either the s-naïve, s-standard or s-modeled paleodistributions at location j. This index is a modified Sørenson

similarity index used to compare two maps of predicted probabilities that each species is present in a given cell rather than relying on presence/absence information. Scores of 1 indicate perfect agreement between the two maps, while an ESP value of 0 indicates complete geographic separation (Godsoe, 2013). We incorporated the uncertainty of each model by randomly permuting the value of each cell 100 times with a Gaussian distribution defined by the mean and standard deviation obtained from model calibration. We report the mean ESP value for each of the s-naïve, s-standard or smodeled paleo-distributions when compared to the hindcast paleo-distribution, as well as ESP values among the three methods. We use repeated measures with species level random-effects to test if the permutated ESP scores were different between the s-naïve, s-standard or s-modeled paleo-distributions. We report marginal F-tests for significance. 2.7. Relationship of sampling bias (s) to species distributions (p) We hypothesized that the degree to which our s-modeled bias correction method could improve estimates of paleo-distributions would be dependent on the congruence between the species distribution (p) and our estimate of sampling bias (s). In species where s and p are similar, there may be less need to precisely quantify s because their combined distribution may sufficiently represent the species distribution (p). We expected that species with a p that was similar to s would show limited improvement from a s-modeled approach and therefore we measured the degree of correlation between each the hindcast distribution (p) of each species and our model of s with Pearson's correlation coefficient. We plot this value (correlation between p and s) against the % improvement of the s-modeled paleo-distribution over the s-naive estimate and use OLS regression to identify any pattern between the two. We hypothesized that species with high correlation between p and s would show reduced improvement from the smodeled bias correction method. 2.8. Validation with pollen We also compared the four paleo-distributions of each species to pollen records from lake sediment cores taken in within 100 km of our study area in the western conterminous USA obtained from Pangea (https://www.pangaea.de) in January 2018 (Appendix G). We selected pollen count data for the genus of each of our focal species with calibrated 14C dates spanning the mid-Holocene. For each sediment core location, we calculated the mean habitat potential value within a 100 km radius for each of the four paleodistributions to account for short-distance wind transport (Broadhurst, 2015). We hypothesized that habitat potential from the s-modeled bias correction method would be higher in the areas surrounding each sediment core containing pollen (dated to the mid-Holocene) than habitat potential from the s-naïve or s-standard methods, and would be more similar to the hindcast distributions. 2.9. Geographic distribution shifts and niche breadth We assumed that dominant niche characteristics have remained stable since the mid-Holocene, and that any changes due to climate would result in distribution shifts only. We therefore investigated how sampling bias could affect estimates of geographic range shifts by measuring changes in spatial predictions of habitat potential between the present-day and mid-Holocene derived from the snaïve, s-standard, s-modeled paleo-distributions. We identified areas where habitat potential either increased or decreased significantly as areas where the difference was greater than twice

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125

environmental conditions) possible in our study area; namely a species occurring on a single grid cell with the most unique environmental conditions. We use the same permutation test described above to estimate 100 realizations of niche breadth for each of the hindcast and paleo-SDM distributions, and use mixed-models with repeated measures to assess differences. We report marginal F-tests for significance.

the combined standard error from the two time periods to capture ~95% of the potential variability due to model error. We report the total area where habitat potential either significantly increased or decreased for each species as estimates of range shifts. We hypothesized that the estimated net change for each species would be most similar between the hindcast and s-modeled paleodistributions. We also posited that not all significant changes in habitat potential would be substantial. For example, areas with habitat potential of 0.95 during the mid-Holocene that changed to 0.85 in the present-day should not constitute substantial changes because these two values are high, even if the total change (0.1) was greater than two times the standard error of each model. We therefore identify areas with substantial change as those where habitat potential was significant and switched from being below 0.5 to above 0.5 (substantial increase), or from being above 0.5 to being below 0.5 (substantial decrease). To test how estimates of environmental niche breadth were affected by sampling bias, we compared estimates of niche breadth for the mid-Holocene for each of the four paleo-distributions. A vast array of metrics for measuring niche breadth have been proposed, and can range from volumetric measurements of n-dimensional hypervolumes (e.g. Blonder et al., 2014), to reduced principle component axes of presumed occupied regions of environmental space (e.g. Saupe et al., 2015), to counts of unique habitat types presumed to be occupied (e.g. Nürnberg and Aberhan, 2013). We developed an approach to quantify the environmental uniqueness of occupied habitat that draws on Mahalanobis distances in environmental space because we were most interested in how unique the occupied habitat was from the rest of the study area across all environmental explanatory variables. The niche breadth value for each species was defined as the median of the squared Mahalanobis distance of all occupied cells:

f2 ¼ ðХ  m Þ0 S1 ðХ  m Þ D o a a a

3. Results 3.1. s-modeled components The hindcast of present-day Neotoma distributions used to develop the availability component of the s-modeled method yielded an AUC score of 0.639 on the 202 FAUNMAP macrofossil test locations, whereas the availability, preservation and discovery components of our s-modeled method resulted in test AUC scores of 0.931, 0.950, and 0.891, respectively. Each component is mapped in Appendix H. 3.2. Present-day and mid-holocene distributions The geographic means of habitat potential scores for presentday conditions in the study area ranged from 0.173 (sd ¼ 0.226; C. ramosissima) to 0.543 (sd ¼ 0.163; A. tridentata), while hindcast paleo-distribution scores ranged from 0.159 (sd ¼ 0.242; C. ramosissima) to 0.552 (sd ¼ 0.157; A. tridentata; Table 1), indicating that habitat potential for each species remained nearly consistent between periods. Mapped predictions for current and mid-Holocene periods are shown in Appendix I. When validated with macrofossils, hindcast paleo-distributions showed test AUC scores ranging from 0.622 (P. ponderosa) to 0.958 (C. ramosissima; Table 1). Of the three bias-correction methods (s-naïve, s-standard and s-modeled) used to calibrate paleo-distributions from macrofossils, the s-naïve resulted in the lowest habitat potential scores across the study region and were substantially lower than hindcast paleo-distributions for all species except Q. gambelii (Table 1). ESP overlap scores between hindcast paleo-distributions and the three fossil calibrated methods (s-naïve, s-standard and s-modeled) were significantly different from one another (F2,1792 ¼ 1274.8605, p < 0.0001), with s-naïve paleo-distributions showing the lowest ESP scores among all species except C. ramosissima (Fig. 1), thereby indicating that without bias correction, sampling bias caused paleo-distributions to be far from their (assumed) true hindcast values. Similarly, the s-modeled paleo-distributions yielded ESP overlap scores that were substantially higher than either the snaïve or s-standard paleo-distributions for all species except Coleogyne ramosissima (Fig. 1), suggesting that the model of s was able to overcome some of the bias inherent in the packrat midden fossil record and bring paleo-distributions closer to their (assumed) true hindcast value. Improvements over the s-standard method

(2)

where Xo is the matrix of explanatory variables used to define the species' niche over all occupied cells, Xa is the matrix of the same explanatory variables over all cells in the study area, ma is a vector of variable means of the study area and S is the covariance matrix of Xa. We determined occupied cells probabilistically, with the equation:

Pi ¼

1

(3)

xi 0:65

1 þ e 0:05

119

where xi is the habitat potential for cell i. Niche breadth values in our study area can range from near 4 (the most generalist species, or organism that has an ability to live in a wide variety of habitats in a wide range of environmental conditions) to over 400, a completely unrealistic value representing the most specialist of species (an organism capable of tolerating a narrow range of

Table 1 Mean (and standard deviation) of habitat potential values in the study area for Present-day (1950e2000 yr) and mid-Holocene (8.4 ka e 4.2 ka) conditions. Habitat potential for the mid-Holocene conditions was estimated using 4 methods: hindcasting from present day conditions (Hindcast), paleo-SDM without bias correction (s-naïve), paleo-SDM with standard correction (s-standard), and paleo-SDM with model correction (s-modeled). Area Under the receiver/operator Curve (AUC) shown in italics for hindcast paleodistributions. Species

Present-Day

Abies concolor Artemisia tridentata Coleogyne ramosissima Juniperus communis Pinus ponderosa Quercus gambelii

0.304 0.543 0.173 0.232 0.332 0.318

mid-Holocene Hindcast (0.224) (0.163) (0.226) (0.265) (0.256) (0.278)

0.273 0.552 0.159 0.245 0.355 0.283

(0.229) (0.157) (0.242) (0.275) (0.268) (0.259)

AUC

snaïve

0.744 0.810 0.958 0.798 0.622 0.758

0.103 0.133 0.210 0.226 0.324 0.457

sstandard (0.146) (0.193) (0.189) (0.196) (0.213) (0.196)

0.120 0.208 0.226 0.258 0.375 0.506

smodeled (0.154) (0.221) (0.188) (0.205) (0.211) (0.150)

0.217 0.378 0.546 0.355 0.541 0.738

(0.299) (0.287) (0.286) (0.318) (0.352) (0.230)

120

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125

Fig. 1. Mean Expected fraction of Shared Presences (ESP) scores for overlap between hindcasting from present day conditions to paleo-SDM without bias correction (s-naïve, light grey), paleo-SDM with standard correction (s-standard, dark grey), and paleo-SDM with model correction (s-modeled, black). ESP measures overlap between s-modeled paleodistribution and hindcast paleo-distribution.

ranged from 24 to 64% with the s-modeled bias correction method (Fig. 1), and up to 93% over s-naïve paleo-distributions. C. ramosissima was the only species that showed a decrease in mean overlap score with the s-modeled bias correction method (Fig. 1); we provide the mean overlap between hindcast and the three paleo-SDM methods (s-naïve, s-standard and s-modeled) in Appendix J. 3.3. Relationship of sampling bias (s) to species distributions (p) Contrary to our hypotheses, we found no evidence that species exhibiting high correlation between p and the estimate of s showed lower improvements (s-modeled over s-naïve) in the ability to recreate hindcast distributions (F1,4 ¼ 2.7562, p ¼ 0.1722). Similarly, we did not find any relationship between macrofossil sample size and improvements in the ability to recreate hindcast distributions (F1,4 ¼ 2.5775, p ¼ 0.1837). In both cases however, C. ramosissima was an anomaly; the s-modeled bias correction method resulted in poor estimates of the mid-Holocene distribution, but also had the highest correlation between p and the estimate of s, as well as the lowest macrofossil sample size (Table 2, Table 2 Percent improvement in the Expected fraction of Shared Presences overlap metric (ESP) of the paleo-SDM with model correction (s-modeled) over the paleo-SDM without bias correction (s-naïve) and paleo-SDM with standard correction (sstandard). Pearson's correlation between the estimate of s and the hindcast (p) paleo-distribution for each species. ESP measures the degree of overlap between two geographic distributions. % Improvement in ESP score Species

s-naïve

s-standard

Abies concolor Artemisia tridentata Coleogyne ramosissima Juniperus communis Pinus ponderosa Quercus gambelii

86.1% 93.1% 26.1% 52.0% 73.4% 31.0%

63.6% 42.1% 24.2% 43.2% 59.9% 24.3%

Correlation between s and p

0.183 0.169 0.530 0.310 0.239 0.053

Appendix B). 3.4. Validation with pollen Independent pollen locations yielded insight consistent with the improvement from the s-modeled method, with five taxa showing higher habitat potential estimated with the s-modeled method at sites where the pollen record indicated presence of the genera during the mid-Holocene (Fig. 2). The sixth species, C. ramosissima, did not have any corresponding pollen records and was excluded. The s-standard method also predicted increased habitat potential (over the s-naïve method) at these sites, but increases were not as great as with the s-modeled method (Fig. 2). 3.5. Geographic distribution shifts and niche breadth Estimates of shifts between the mid-Holocene and present-day varied dramatically among the methods of estimating paleodistributions. Projections of areas with significant increases of habitat potential (range expansion) from the mid-Holocene to present-day with s-modeled paleo-distributions provided the closest estimates to changes in habitat potential from hindcast distributions for all species except A. concolor and Q. gambelii ( Table 3), where the s-standard paleo-distribution provided the better estimate. However, when assessing significant reductions in habitat potential (range decreases), the s-modeled paleo-distributions did not improve estimates over the s-naïve paleo-distributions (Table 4). This suggests that while the s-modeled paleo-distributions provided higher spatial congruence with hindcast distributions than the snaïve or s-standard paleo-distributions, they may overestimate habitat potential in the mid-Holocene and therefore result in increased estimates of habitat loss from the mid-Holocene to the present. Hindcast niche breadths for the 6 species ranged from 5.59 to 10.02 (Table 5), with A. concolor and J. communis being the most specialized, and A. tridentata and Q. gambelii the most generalist species within our study area. Niche breadths estimated with the s-

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125

121

Fig. 2. Mean habitat potential values within 100 km radius of sites where the pollen record indicated presence of the genera during the mid-Holocene for three paleo-SDM distributions: paleo-SDM without bias correction (s-naïve), paleo-SDM with standard correction (s-standard), and paleo-SDM with model correction (s-modeled).

modeled paleo-distributions most closely replicated those from hindcast distributions for three species; the s-standard method best replicated hindcast niche breadths for one species, and the snaïve best replicated hindcast niche breadth for two (Table 5). The least accurate estimate of hindcast niche breadth came from the smodeled paleo-distribution, which was 26.1% from the assumed true value for Q. gambelii. Estimates of niche breadths for the three paleo-SDM methods were significantly different from one anther (F2,352 ¼ 56.8256, p ¼ 0.0196) when accounting for differences among species, but the s-standard and s-modeled estimates were not (F1,233 ¼ 1.3003, p ¼ 0.2553). 4. Discussion 4.1. Packrat midden macrofossil bias Paleobiogeographers have studied fossilized packrat middens since the early work of Philip Wells in the 1960s (e.g. Wells, 1966). In the last half-century, numerous researchers have contributed to packrat midden archives, which have continued to grow taxonomically, temporally and spatially (Williams et al., 2018). These archives offer access to a wealth of macrofossil observations, many of which have species level taxonomic precision and can be geo-

located to within 100 m of the midden site from which they were extracted (Vaughan, 1990). Observations span well into the Pleistocene at the limits of 14C dating methods, making them ideally suited for investigating changes in recent north American floral (e.g. Cole and Webb, 1985; Jackson et al., 2005; McAuliffe and Van Devender, 1998; Spaulding, 1990; Thompson and Anderson, 2000), faunal (e.g. Mead, 1981; Van Devender and Mead, 1978; Van Devender et al., 1977) and climatic (e.g. Arundel, 2002; Coats et al., 2008; Jacobson, 1991; Thompson et al., 2012) reconstructions. While the breadth of biogeographic inquiries using these data continues to grow (Williams et al., 2018), their use for paleo-SDM has often been limited, likely due to the spatial bias inherent in the NAPM database (Mensing et al., 2000; Webb and Betancourt, 1990). In contrast, pollen from sediment cores and other fossil distribution data have been used widely for paleo-SDM (Maguire et al., 2015; Moreno-Amat et al., 2017; Myers et al., 2015; Svenning et al., 2011; Varela et al., 2011) and have led to investigations of niche stability (e.g. Saupe et al., 2014; Stigall, 2012; Veloz et al., 2012), glacial & interglacial extinctions (e.g. Lorenzen et al., 2011), mammalian community structure (e.g. Rowan et al., s-Bravo et al., 2008; Veloz 2016), range dynamics (e.g. Nogue ri, 2008), and conet al., 2012), speciation (e.g. Peterson and Nya servation priorities (e.g. Williams et al., 2013), among others.

Table 3 Total area (km2) with significant increase in habitat potential between mid-Holocene and present-day conditions. Areas with a significant increase in habitat potential were defined as areas where the difference was greater than twice the combined standard error of the paleo-distributions from the two time periods. Areas with substantial increase in habitat potential were identified as areas that were significant and where habitat potential switched from being below 0.5 to above 0.5 (substantial increase). Study area was 3,162,474 km2 in size. species

Range Expansion (Increase) significant

Abies concolor Artemisia tridentata Coleogyne ramosissima Juniperus communis Pinus ponderosa Quercus gambelii

substantial

Hindcast

s-naïve

s-standard

s-modeled

Hindcast

s-naïve

s-standard

s-modeled

360,035 10,819 240,286 11,933 11,699 109,392

1,764,634 2,602,924 65,178 328,123 612,110 292,104

1,539,659 2,302,023 52,025 291,404 522,412 155,081

2,088,242 2,139,835 281,451 90,311 266,302 19,230

79,818 456 28,653 4091 4257 33,186

487,001 1,493,239 42,955 223,126 523,723 261,796

457,162 1,239,031 40,020 203,846 469,307 109,923

430,001 1,113,236 155,056 52,362 48,367 168

122

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125

Table 4 Total area (km2) with significant decrease in habitat potential between mid-Holocene and present-day conditions. Areas with a significant decrease in habitat potential were defined as areas where the difference was less than the negative of twice the combined standard error of the paleo-distributions from the two time periods. Areas with substantial decrease in habitat potential were identified as areas that were significant where habitat potential switched from being above 0.5 to below 0.5 (substantial decrease). species

Range Contraction (Decrease) significant

Abies concolor Artemisia tridentata Coleogyne ramosissima Juniperus communis Pinus ponderosa Quercus gambelii

substantial

Hindcast

s-naïve

s-standard

s-modeled

Hindcast

s-naïve

s-standard

s-modeled

45,006 85,029 59,402 105,601 227,459 10,189

18,767 18,936 42,186 118,229 397,887 1,089,216

20,329 24,692 50,747 256,707 990,500 1,174,115

502,363 506,071 2,424,022 2,146,203 2,168,690 2,849,691

5828 19,144 21,048 9453 93,143 3664

16,101 586 29,028 74,457 205,150 741,152

17,300 1904 32,316 93,864 310,832 847,762

2,849,691 76,783 1,466,654 264,615 807,995 1,712,243

Table 5 Niche breadth as estimated by each of the 4 methods: hindcasting from present day conditions (Hindcast), paleo-SDM without bias correction (s-naïve), paleo-SDM with standard correction (s-standard), and paleo-SDM with model correction (s-modeled). Niche breadth used here measures the uniqueness of the environmental conditions defining each species' geographic distribution and was defined as the median of the squared Mahalanobis distance of all occupied cells. Absolute percent difference from hindcast niche breadth shown in parentheses. species

Abies concolor Artemisia tridentata Coleogyne ramosissima Juniperus communis Pinus ponderosa Quercus gambelii

s-naïve

Hindcast

s-standard

s-modeled

mean

mean

% Diff

mean

% Diff

mean

% Diff

10.02 5.59 6.27 9.52 7.74 5.74

12.05 4.56 7.86 8.61 6.59 4.80

(20.2) (18.3) (25.3) (9.5) (14.9) (16.3)

11.25 4.63 7.53 8.55 6.47 4.80

(12.2) (17.1) (20.2) (10.2) (16.5) (16.5)

10.37 6.99 6.92 8.35 6.83 4.25

(3.5) (24.9) (10.4) (12.3) (11.7) (26.1)

We found that paleo-distributions derived from hindcasting provided reasonably accurate reconstructions of mid-Holocene distributions when assessed with macrofossils obtained from the NAPM database. We therefore used these hindcast distributions as our reference paleo-distributions from which we evaluated multiple methods of bias correction used in paleo-SDM. We found that paleo-distributions for the mid-Holocene modeled without bias correction (s-naïve) provided poor matches to hindcast paleodistributions. The s-modeled correction method improved paleodistributions for our six species and showed, on average, 52% improvement in overlap with hindcast paleo-distributions than snaïve paleo-distributions. These improvements were confirmed at independent locations from lake sediment pollen records, where we found habitat potential scores for s-modeled paleodistributions to be higher than the s-naïve and s-standard paleodistributions. We used hindcast paleo-distributions as our reference distributions because they draw on a wealth of biogeographic repositories characterizing contemporary distributions of extant species. Hindcast distributions however, assume that a species niche has remained unchanged between the calibration and projection time periods. This assumption is reasonable for periods in the mid or late Holocene, an interval of only a few millennia when the world's climate was relatively stable and substantial evolutionary change unlikely, but may be less so for paleobiogeographic studies of the Pleistocene or earlier (Peterson, 2011; Walls and Stigall, 2011). In studies where hindcasting is less appropriate, paleo-SDM calibrated with paleo-ecological archives has been shown to be useful for estimating paleo-distributions (e.g. Walls and Stigall, 2012). 4.2. Niche characteristics Our work confirms the difficulty of measuring niche breadth by showing how three different methods of paleo-SDM, all calibrated with the same macrofossil records can result in widely varying

estimates of niche breadth. Further, our work suggests that fossil bias may cause distortions in assessments of niche breadth and niche stability through time. Niche stability assumes that fundamental niche characteristics, such as niche breadth, are conserved through time (Wiens and Graham, 2005) and that any changes in climate or other environmental determinants of habitat may instead cause shifts in geographic distributions. Formal tests, such as niche equivalency (Warren et al., 2008), have been developed and used extensively to investigate niche stability in related extant taxa (Saupe et al., 2014; e.g. Strubbe et al., 2014), but these tests compare spatial predictions of habitat potential to infer similarity in niche characteristics, and do not estimate niche properties directly. Metrics of niche breadth can span volumetric measurements of n-dimensional hypervolumes or dominant principle component axes of presumed occupied regions of environmental space to counts of unique habitat types presumed to be occupied. These different niche breadth metrics, including ours, are highly dependent on the environmental factors used to define an organism's niche, and as such, are often not comparable (Slatyer et al., 2013). We estimated niche breadth using a common set of environmental factors to make estimates comparable among species, but found that the s-modeled bias correction method improved niche breadth estimates in only three of the six species. 4.3. Taxonomic bias While our work focuses on spatial bias inherent in the NAPM database, taxonomic bias is another form often faced when leveraging paleoecological archives. In its most basic form, taxonomic bias can result in a greater prevalence of fossils for some species, but not others; this can lead to biased inference of abundance and diversity. This is especially true when incorporating macrofossils obtained from packrat middens where vegetative material from species such as Fallugia paradoxa, Juniperus monosperma, Yucca baileyi, and Ephedra torreyana, among others, may be

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125

collected at higher frequencies than surrounding environs for forage (Dial and Czaplewski, 1990). Other species may be less well represented due to changes in packrat foraging behavior proximal to a midden site (Spaulding et al., 1990) or differences among packrat species in diet preferences (Vaughan, 1990). These differences certainly influence estimates of abundance and diversity, and have been discussed at length elsewhere (c.f. Betancourt et al., 1990). Among our six focal species, at least two (J. communis and P. ponderosa) are likely found in the diets of several packrat species (Vaughan, 1990). These two species also had the highest number of observations in the NAPM database, but did not stand out among our species in terms of relative gains from the s-modeled bias correction over the s-naïve paleo-distributions. Moreover, our goal was not to estimate paleo-abundances or indices of diversity among taxa, but rather, to offer and test an alternative method to account for spatial bias in paleoecological archives when using paleo-SDM.

123

adapted for other paleoecological archives to investigate finer time slices and earlier periods of our earth's history. Data availability All relevant data are hosted by ScienceBase, a public repository hosted by the U.S. Geological Survey. Data can be accessed at the following DOI: https://doi.org/10.5066/P9843JFT. Authors contributions RI and JF conceived the ideas and designed methodology; RI collected and analyzed the data; RI led the writing of the manuscript, JF, TE and KN contributed to the writing and editing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication. Declarations

5. Temporal uncertainty Declarations of interest: none. We incorporated estimates of model uncertainty in our permutation tests of ESP scores, but were unable to include measures of temporal uncertainty from the calibrated ages of macrofossils or the downscaled paleoclimate simulations. Instead, we limited our analyses to the Northgrippian subdivision of the Holocene (8.3 ka e 4.2 ka), an interval of only 4 ka centered around the 6 ka experimental design for mid-Holocene boundary conditions used to simulate mid-Holocene climate (Braconnot et al., 2012). This interval, while representative of the downscaled paleoclimate simulations, dramatically constrains available sample sizes from the NAPM database; any further restrictions to the study period would limit species available for consideration. Our smallest sample size for model calibration was 4, the bare minimum needed to calibrate distribution models for narrow niched species (van Proosdij et al., 2015). Paleoecological archives, such as the European Pollen Database (http://www.europeanpollendatabase.net), may offer far greater sample sizes, thereby allowing narrow time intervals for selecting calibration records (Moreno-Amat et al., 2015). 6. Access to paleoecological archives The marked increase in online data storage and ease of establishing and maintaining web services has led to the proliferation of online, open access paleoecological databases well suited for paleoSDM. Such resources include broad data repositories like PANGAEA (https://pangaea.de), and the NOAA/World Data Service for Paleoclimatology (https://www.ncdc.noaa.gov), as well as curated databases such as the NAPM database (https://geochange.er.usgs.gov/ midden/; Strickland et al., 2013), the Neotoma Paleocology Database (Williams et al., 2018), FAUNMAP and MIOMAP (Carrasco et al., 2005), among others. These databases and the software to access them (e.g. Varela et al., 2014) offer spatially referenced fossil archives that extend well into the Miocene. When species paleodistribution records are combined with modeled and reconstructed spatial paleoclimate simulations such as the Coupled Modeling Intercomparison Project (CMIP5: http://cmip-pcmdi.llnl. gov/cmip5/) and the Paleoclimate Modeling Intercomparison Project (PMIP3: https://pmip3.lsce.ipsl.fr/) (e.g. Fordham et al., 2017; Lima-Ribeiro et al., 2015; Lorenz et al., 2016), paleo-SDM offers biogeographers new tools for exploring paleo-distributions well beyond the mid-Holocene addressed here. Moreover, improvements in the temporal resolution of gridded climate data offer opportunities to investigate many time slices of recent history (Blois et al., 2013; Pearman et al., 2008). We suggest that the approach described here to estimate sampling bias (s) can be

Acknowledgments Funding for this work was provided by the USGS Ecosystems Mission Area and Western Ecological Research Center. We thank R. Dorn and P. Henne for reviewing initial drafts of this work and L. Strickland for providing access and assistance with the NAPM database. We also thank 2 anonymous reviewers for their valuble comments on this manuscript. Any use of trade, product or firm names in this publication is for descriptive purposes only and does not imply endorsement by the US goverment. Appendix A. Supplementary data Supplementary data related to this article can be found at https://doi.org/10.1016/j.quascirev.2018.08.015. References Allison, P.A., Bottjer, D.J., 2010. Taphonomy. Springer. Angulo, D.F., Amarilla, L.D., Anton, A.M., Sosa, V., 2017. Colonization in North american arid lands: the journey of agarito (berberis trifoliolata) revealed by multilocus molecular data and packrat midden fossil remains. PLoS One 12. https://doi.org/10.1371/journal.pone.0168933 e0168933e25. Araújo, M.B., Guisan, A., 2006. Five (or so) challenges for species distribution modelling. J. Biogeogr. 33, 1677e1688. https://doi.org/10.1111/j.13652699.2006.01584.x. Arundel, S.T., 2002. Modeling climate limits of plants found in sonoran desert packrat middens. Quat. Res. (Tokyo) 58, 112e121. https://doi.org/10.1006/ qres.2002.2349. Baddeley, A., Rubak, E., Turner, R., 2015. Spatial Point Patterns: Methodology and Applications with R. Chapman and Hall/CRC Press, London. Betancourt, J.L., Van Devender, T.R., Martin, P.S., 1990. Packrat Middens: the Last 40,000 Years of Biotic Change. University of Arizona Press, Tucson. Blois, J.L., Williams, J.W., Fitzpatrick, M.C., Jackson, S.T., Ferrier, S., 2013. Space can substitute for time in predicting climate-change effects on biodiversity. Proc. Natl. Acad. Sci. U. S. A. 110, 9374e9379. https://doi.org/10.1073/ pnas.1220228110. Blonder, B., Lamanna, C., Violle, C., Enquist, B.J., 2014. The n-dimensional hypervolume. Global Ecol. Biogeogr. 23, 595e609. https://doi.org/10.1111/geb.12146. Boria, R.A., Olson, L.E., Goodman, S.M., Anderson, R.P., 2014. Spatial filtering to reduce sampling bias can improve the performance of ecological niche models. Ecol. Model. 275, 73e77. https://doi.org/10.1016/j.ecolmodel.2013.12.012. Braconnot, P., Harrison, S.P., Kageyama, M., Bartlein, P.J., Masson-Delmotte, V., AbeOuchi, A., Otto-Bliesner, B., Zhao, Y., 2012. Evaluation of climate models using palaeoclimatic data. Nat. Clim. Change 2, 417e424. https://doi.org/10.1038/ nclimate1456. Brewer, S., Jackson, S.T., Williams, J.W., 2012. Paleoecoinformatics: applying geohistorical data to ecological questions. Trends Ecol. Evol. 27, 104e112. https:// doi.org/10.1016/j.tree.2011.09.009. Broadhurst, L., 2015. Pollen dispersal in fragmented populations of the dioecious wind-pollinated tree, Allocasuarina verticillata (drooping sheoak, drooping sheoak; allocasuarinaceae). PLoS One 10. https://doi.org/10.1371/

124

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125

journal.pone.0119498 e0119498e17. Brotons, L., Thuiller, W., Araujo, M., Hirzel, A., 2004. Presence-absence versus presence-only modelling methods for predicting bird habitat suitability. Ecography 27, 437e448. Carnaval, A.C., Moritz, C., 2008. Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest. J. Biogeogr. 35, 1187e1201. https://doi.org/10.1111/j.1365-2699.2007.01870.x. Carrasco, M.A., Kraatz, E.B., Davis, E.B., Barnosky, A.D., 2005. Miocene Mammal Mapping Project (MIOMAP). Coats, L.L., Cole, K.L., Mead, J.I., 2008. 50,000 years of vegetation and climate history on the Colorado Plateau, Utah and Arizona, USA. Quat. Res. (Tokyo) 70, 322e338. https://doi.org/10.1016/j.yqres.2008.04.006. Cohen, K.M., Finney, S.C., Gibbard, P.L., Fan, J.X., 2013. The ICS International Chronostratigraphic Chart, p. 199 (Episodes). Cole, K.L., Webb, R.H., 1985. Late Holocene vegetation changes in greenwater valley, mojave desert, California. Quat. Res. (Tokyo) 23, 227e235. https://doi.org/ 10.1016/0033-5894(85)90030-4. Danielson, J.J., Gesch, D.B., 2011. Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010). U.S. Geological Survey Open-file. Report 2011e1073. Dial, K.P., Czaplewski, N.J., 1990. Do woodrat middens accurately represent the animals' environments and diets? The Woodhouse Mesa study. In: Packrat Middens: the Last 40,000 Years of Biotic Change. University of Arizona Press, Tucson. Elith, J., Phillips, S.J., Hastie, T., Dudik, M., Chee, Y.E., Yates, C.J., 2011. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 17, 43e57. https://doi.org/ 10.1111/j.1472-4642.2010.00725.x. Fielding, A., Bell, J., 1997. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 24, 38e49. , F., Haythorne, S., Wigley, T.M.L., Otto-Bliesner, B.L., Chan, K.C., Fordham, D.A., Saltre Brook, B.W., 2017. PaleoView: a tool for generating continuous climate projections spanning the last 21,000 years at regional and global scales. Ecography 1e29. https://doi.org/10.1111/ecog.03031. Franklin, J., 2010. Mapping Species Distributions. Cambridge University Press. Franklin, J., Potts, A.J., Fisher, E.C., Cowling, R.M., Marean, C.W., 2015. Paleodistribution modeling in archaeology and paleoanthropology. Quat. Sci. Rev. 110, 1e14. https://doi.org/10.1016/j.quascirev.2014.12.015. Gent, P.R., Danabasoglu, G., Donner, L.J., Holland, M.M., Hunke, E.C., Jayne, S.R., Lawrence, D.M., Neale, R.B., Rasch, P.J., Vertenstein, M., Worley, P.H., Yang, Z.-L., Zhang, M., 2011. The community climate System model version 4. J. Clim. 24, 4973e4991. https://doi.org/10.1175/2011JCLI4083.1. Gill, J.L., Williams, J.W., Jackson, S.T., Lininger, K.B., Robinson, G.S., 2009. Pleistocene megafaunal collapse, novel plant communities, and enhanced fire regimes in North America. Science 326, 1100e1103. https://doi.org/10.1126/ science.1179504. Godsoe, W., 2013. Inferring the similarity of species distributions using Species' Distribution Models. Ecography 37, 130e136. https://doi.org/10.1111/j.16000587.2013.00403.x. Haslett, J., Parnell, A., 2008. A simple monotone process with application to radiocarbon-dated depth chronologies. J Roy Stat Soc C-App 57, 399e418. https://doi.org/10.1111/j.1467-9876.2008.00623.x. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965e1978. https://doi.org/10.1002/joc.1276. Jackson, S.T., Betancourt, J.L., Lyford, M.E., Gray, S.T., Rylander, K.A., 2005. A 40,000year woodrat-midden record of vegetational and biogeographical dynamics in north-eastern Utah, USA. J. Biogeogr. 32, 1085e1106. https://doi.org/10.1111/ j.1365-2699.2005.01251.x. Jacobson, G.L., 1991. Packrat middens and climate change. Ecology 72. https:// doi.org/10.2307/2937220, 760e760. lez-Hern Lima-Ribeiro, M.S., Varela, S., Gonza andez, J., de Oliveira, G., DinizFilho, J.A.F., Terribile, L.C., 2015. EcoClimate: a database of climate data from multiple models for past, present, and future for macroecologists and biogeographers. Biodivers. Inf. 10 https://doi.org/10.17161/bi.v10i0.4955. Lorenz, D.J., Nieto-Lugilde, D., Blois, J.L., Fitzpatrick, M.C., Williams, J.W., 2016. Downscaled and debiased climate simulations for North America from 21,000 years ago to 2100AD. Sci. Data 3. https://doi.org/10.1038/sdata.2016.48, 160048e19. s-Bravo, D., Orlando, L., Weinstock, J., Binladen, J., Lorenzen, E.D., Nogue Marske, K.A., Ugan, A., Borregaard, M.K., Gilbert, M.T.P., Nielsen, R., Ho, S.Y.W., Goebel, T., Graf, K.E., Byers, D., Stenderup, J.T., Rasmussen, M., Campos, P.F., Leonard, J.A., Koepfli, K.-P., Froese, D., Zazula, G., Stafford, T.W., AarisSørensen, K., Batra, P., Haywood, A.M., Singarayer, J.S., Valdes, P.J., Boeskorov, G., Burns, J.A., Davydov, S.P., Haile, J., Jenkins, D.L., Kosintsev, P., Kuznetsova, T., Lai, X., Martin, L.D., McDonald, H.G., Mol, D., Meldgaard, M., Munch, K., Stephan, E., Sablin, M., Sommer, R.S., Sipko, T., Scott, E., Suchard, M.A., Tikhonov, A., Willerslev, R., Wayne, R.K., Cooper, A., Hofreiter, M., Sher, A., Shapiro, B., Rahbek, C., Willerslev, E., 2011. Species-specific responses of Late Quaternary megafauna to climate and humans. Nature 479, 359e364. https:// doi.org/10.1038/nature10574. Maguire, K.C., Nieto-Lugilde, D., Fitzpatrick, M.C., Williams, J.W., Blois, J.L., 2015. Modeling species and community responses to past, present, and future episodes of climatic and ecological change. Annu. Rev. Ecol. Evol. Syst. 46, 343e368. https://doi.org/10.1146/annurev-ecolsys-112414-054441. Martinez-Meyer, E., Townsend Peterson, A., Hargrove, W.W., 2004. Ecological niches as stable distributional constraints on mammal species, with implications for

Pleistocene extinctions and climate change projections for biodiversity. Global Ecol. Biogeogr. 13, 305e314. https://doi.org/10.1111/j.1466-822X.2004.00107.x. McAuliffe, J.R., Van Devender, T.R., 1998. A 22,000-year record of vegetation change in the north-central Sonoran Desert. Palaeogeogr. Palaeoclimatol. Palaeoecol. 141, 253e275. https://doi.org/10.1016/S0031-0182(98)00054-6. Mead, J.I., 1981. The last 30,000 Years of faunal history within the grand-canyon, Arizona. Quat. Res. (Tokyo) 15, 311e326. https://doi.org/10.1016/00335894(81)90033-8. Mensing, S.A., Elston Jr., R.G., Raines, G.L., Tausch, R.J., 2000. A GIS model to predict the location of fossil packrat (Neotoma) middens in central Nevada. Western North American. https://doi.org/10.2307/41717022. Moreno-Amat, E., Mateo, R.G., Nieto-Lugilde, D., Morueta-Holme, N., Svenning, J.-C., García-Amorena, I., 2015. Impact of model complexity on cross-temporal transferability in Maxent species distribution models: an assessment using paleobotanical data. Ecol. Model. 312, 308e317. https://doi.org/10.1016/ j.ecolmodel.2015.05.035. Moreno-Amat, E., Rubiales, J.M., Morales-Molino, C., García-Amorena, I., 2017. Incorporating plant fossil data into species distribution models is not straightforward: pitfalls and possible solutions. Quat. Sci. Rev. 170, 56e68. https://doi.org/10.1016/j.quascirev.2017.06.022. Myers, C.E., Stigall, A.L., Lieberman, B.S., 2015. PaleoENM: applying ecological niche modeling to the fossil record. Paleobiology 41, 226e244. https://doi.org/ 10.1017/pab.2014.19. s-Bravo, D., 2009. Predicting the past distribution of species climatic niches. Nogue Global Ecol. Biogeogr. 18, 521e531. https://doi.org/10.1111/j.14668238.2009.00476.x. s-Bravo, D., Rodríguez, J., Hortal, J., Batra, P., Araújo, M.B., 2008. Climate Nogue change, humans, and the extinction of the woolly mammoth. PLoS Biol. 6, e79. https://doi.org/10.1371/journal.pbio.0060079. Nürnberg, S., Aberhan, M., 2013. Habitat breadth and geographic range predict diversity dynamics in marine Mesozoic bivalves. Paleobiology 39, 360e372. https://doi.org/10.1666/12047. Ordonez, A., 2013. Realized climatic niche of North American plant taxa lagged behind climate during the end of the Pleistocene. Am. J. Bot. 100, 1255e1265. https://doi.org/10.3732/ajb.1300043. Pearman, P.B., Randin, C.F., Broennimann, O., Vittoz, P., van der Knaap, W.O., Engler, R., Le Lay, G., Zimmermann, N.E., Guisan, A., 2008. Prediction of plant species distributions across six millennia. Ecol. Lett. 11, 357e369. https:// doi.org/10.1111/j.1461-0248.2007.01150.x. Peterson, A.T., 2011. Ecological niche conservatism: a time-structured review of evidence. J. Biogeogr. 38, 817e827. https://doi.org/10.1111/j.13652699.2010.02456.x.  ri, A.S., Peterson, A.T., Nya 2008. Ecological niche conservatism and pleistocene refugia in the thrush-like Mourner, Schiffornis sp., in the neotropics. Evolution 62, 173e183. https://doi.org/10.1111/j.1558-5646.2007.00258.x. Phillips, S.J., Anderson, R.P., Schapire, R.E., 2006. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231e259. https://doi.org/ 10.1016/j.ecolmodel.2005.03.026. Phillips, S.J., Dudik, M., Elith, J., Graham, C., Lehmann, A., Leathwick, J., Ferrier, S., 2009. Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecol. Appl. 19, 181e197. https:// doi.org/10.1890/07-2153.1. R Core Team, 2016. R: a Language and Environment for Statistical Computing. Rowan, J., Kamilar, J.M., Beaudrot, L., Reed, K.E., 2016. Strong influence of palaeoclimate on the structure of modern African mammal communities. Proc. Biol. Sci. 283 https://doi.org/10.1098/rspb.2016.1207, 20161207e9. Saupe, E.E., Hendricks, J.R., Portell, R.W., Dowsett, H.J., Haywood, A., Hunter, S.J., Lieberman, B.S., 2014. Macroevolutionary consequences of profound climate change on niche evolution in marine molluscs over the past three million years. Proc. Biol. Sci. 281 https://doi.org/10.1098/rspb.2014.1995, 20141995e20141995.  n, J., Saupe, E.E., Qiao, H., Hendricks, J.R., Portell, R.W., Hunter, S.J., Sobero Lieberman, B.S., 2015. Niche breadth and geographic range size as determinants of species survival on geological time scales. Global Ecol. Biogeogr. 24, 1159e1169. https://doi.org/10.1111/geb.12333. Schoener, T.W., 1968. The Anolis lizards of bimini: resource partitioning in a complex fauna. Ecology 49, 704e726. https://doi.org/10.2307/1935534. Slatyer, R.A., Hirst, M., Sexton, J.P., 2013. Niche breadth predicts geographical range size: a general ecological pattern. Ecol. Lett. 16, 1104e1114. https://doi.org/ 10.1111/ele.12140. Spaulding, W.G., 1990. Vegetation dynamics during the last deglaciation, southeastern Great Basin, U.S.A. Quat. Res. (Tokyo) 33, 188e203. https://doi.org/ 10.1016/0033-5894(90)90018-g. Spaulding, W.G., Betancourt, J., Croft, L.K., Cole, K.L., 1990. Packrat middens: their composition and methods of analysis. In: Packrat Middens: the Last 40,000 Years of Biotic Change. University of Arizona Press, Tucson. Stigall, A.L., 2012. Using ecological niche modelling to evaluate niche stability in deep time. J. Biogeogr. 39, 772e781. https://doi.org/10.1111/j.13652699.2011.02651.x. Strickland, L.E., Thompson, R.S., Anderson, K.H., 2013. Usgs/noaa North American Packrat Midden Database, Data Dictionary. BiblioGov. Strubbe, D., Beauchard, O., Matthysen, E., 2014. Niche conservatism among nonnative vertebrates in Europe and North America. Ecography 38, 321e329. https://doi.org/10.1111/ecog.00632. s-Bravo, D., Normand, S., 2011. Svenning, J.-C., Fløjgaard, C., Marske, K.A., Nogue

R. Inman et al. / Quaternary Science Reviews 198 (2018) 115e125 Applications of species distribution modeling to paleobiology. Quat. Sci. Rev. 30, 2930e2947. https://doi.org/10.1016/j.quascirev.2011.06.012. Swetnam, T.W., Allen, C.D., 1999. Applied historical ecology: using the past to manage for the future. Ecol. Appl. 9, 1189e1206. https://doi.org/10.1890/10510761(1999)009 [1189:aheutp]2.0.co;2. Sørensen, T., 1948. A Method of Establishing Groups of Equal Amplitude in Plant Sociology Based on Similarity of Species Content and its Application to Analyses of the Vegetation on Danish Commons. Thompson, R.S., Anderson, K.H., 2000. Biomes of western North America at 18,000, 6000 and 0 14C yr bp reconstructed from pollen and packrat midden data. J. Biogeogr. 27, 555e584. https://doi.org/10.1046/j.1365-2699.2000.00427.x. Thompson, R.S., Anderson, K.H., Pelltier, R.T., Strickland, L.E., Bartlein, P.J., Shafer, S.L., 2012. Quantitative estimation of climatic parameters from vegetation data in North America by the mutual climatic range technique. Quat. Sci. Rev. 51, 18e39. https://doi.org/10.1016/j.quascirev.2012.07.003. Van Devender, T.R., Mead, J.I., 1978. Early Holocene and late Pleistocene amphibians and reptiles in Sonoran Desert packrat middens. Copeia 464. https://doi.org/ 10.2307/1443613, 1978. Van Devender, T.R., Phillips, A.M., Mead, J.I., 1977. Late Pleistocene Reptiles and Small Mammals from the Lower Grand Canyon of Arizona. SW. Nat. 22, 49. https://doi.org/10.2307/3670463. van Proosdij, A.S.J., Sosef, M.S.M., Wieringa, J.J., Raes, N., 2015. Minimum required number of specimen records to develop accurate species distribution models. Ecography 39, 542e552. https://doi.org/10.1111/ecog.01509. Varela, S., Gonz alez-Hern andez, J., Sgarbi, L.F., Marshall, C., Uhen, M.D., Peters, S., McClennen, M., 2014. paleobioDB: an R package for downloading, visualizing and processing data from the Paleobiology Database. Ecography 38, 419e425. https://doi.org/10.1111/ecog.01154. Varela, S., Lobo, J.M., Hortal, J., 2011. Using species distribution models in paleobiogeography: A matter of data, predictors and concepts. Palaeogeogr. Palaeoclimatol. Palaeoecol. 310, 451e463. https://doi.org/10.1016/ j.palaeo.2011.07.021. Vaughan, T.A., 1990. Ecology of living packrats. In: Packrat Middens: the Last 40,000 Years of Biotic Change. University of Arizona Press, Tucson. Veloz, S.D., Williams, J.W., Blois, J.L., He, F., Otto-Bliesner, B., Liu, Z., 2012. No-analog climates and shifting realized niches during the late quaternary: implications for 21st-century predictions by species distribution models. Global Change Biol.

125

18, 1698e1713. https://doi.org/10.1111/j.1365-2486.2011.02635.x. Vilhena, D.A., Smith, A.B., 2013. Spatial Bias in the Marine Fossil Record. PLoS One 8, e74470ee74478. https://doi.org/10.1371/journal.pone.0074470. Walls, B.J., Stigall, A.L., 2012. A field-based analysis of the accuracy of niche models applied to the fossil record. Paleontological Contributions. https://doi.org/ 10.17161/PC.1808.8950. Walls, B.J., Stigall, A.L., 2011. Analyzing niche stability and biogeography of Late Ordovician brachiopod species using ecological niche modeling. Palaeogeogr. Palaeoclimatol. Palaeoecol. 299, 15e29. https://doi.org/10.1016/ j.palaeo.2010.10.024. Warren, D.L., Glor, R.E., Turelli, M., 2008. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62, 2868e2883. https://doi.org/10.1111/j.1558-5646.2008.00482.x. Webb, R.H., Betancourt, J.L., 1990. The spatial and temporal distribution of radiocarbon ages from pack-rat middens. In: Packrat Middens: the Last 40,000 Years of Biotic Change. University of Arizona Press, Tucson. Wells, P.V., 1966. Late Pleistocene Vegetation and Degree of Pluvial Climatic Change in the Chihuahuan Desert. Science 153, 970e975. https://doi.org/10.1126/ science.153.3739.970. Wiens, J., Graham, C., 2005. Niche Conservatism: Integrating Evolution, Ecology, and Conservation Biology. Annu. Rev. Ecol. Evol. Syst. 36, 519e539. https://doi.org/ 10.1146/annurev.ecolsys.36.102803.095431. Wiens, J.A., Stralberg, D., Jongsomjit, D., Howell, C.A., Snyder, M.A., 2009. Niches, models, and climate change: assessing the assumptions and uncertainties. Proc. Natl. Acad. Sci. Unit. States Am. 106 (Suppl. 2), 19729e19736. https://doi.org/ 10.1073/pnas.0901639106. Williams, J.W., Grimm, E.C., Blois, J.L., Charles, D.F., Davis, E.B., Goring, S.J., Graham, R.W., Smith, A.J., Anderson, M., Arroyo-Cabrales, J., Ashworth, A.C., Betancourt, J.L., Bills, B.W., Booth, R.K., Buckland, P.I., Curry, B.B., Giesecke, T., Jackson, S.T., Latorre, C., Nichols, J., Purdum, T., Roth, R.E., Stryker, M., Takahara, H., 2018. The Neotoma Paleoecology Database, a multiproxy, international, community-curated data resource. Quat. Res. 89, 156e177. https:// doi.org/10.1017/qua.2017.105. Williams, J.W., Kharouba, H.M., Veloz, S., Vellend, M., McLachlan, J., Liu, Z., OttoBliesner, B., He, F., 2013. The ice age ecologist: testing methods for reserve prioritization during the last global warming. Global Ecol. Biogeogr. 22, 289e301. https://doi.org/10.1111/j.1466-8238.2012.00760.x.