Population reconstructions for humans and megafauna suggest mixed causes for North American Pleistocene extinctions | Nature Communications

Compiling radiocarbon datasets

We compiled 521 direct radiocarbon dates on extinct megafauna derived from published sources for the contiguous US including a small sample (n = 17) of specimens derived from the southern margins of directly adjacent Canadian provinces (i.e., Alberta, Ontario, and Nova Scotia; Fig. 1; Supplementary Data 1). We included only direct AMS or conventional radiocarbon assays on hide, hoof, dung, and bone collagen; we excluded dates on apatite, whole bone, or bone mineral. In this comprehensive megafauna radiocarbon dataset we included genera that went extinct in North America as well as several well-defined extinct species that belong to genera that live on in the region, namely American lion, dire wolf, and Harrington’s mountain goat. So as not to oversample megafauna during the period they overlap with people, we followed Surovell et al.14 and excluded dates derived from unequivocal archeological contexts. To this end, we followed the vetting of Grayson and Meltzer32 and excluded megafaunal dates from sites that represent clear kill and/or scavenging associations.

To obtain archeological radiocarbon dates, we queried the CARD58 for all dates falling within the study region (south of 52 °N latitude) between 15.0 and 10.0 ka, the former chosen as the most widely accepted earliest date for the initial colonization of the Americas (following Weitzel and Codding64). Given that we made use of minimally vetted and user-supplied radiocarbon data from the CARD, it was necessary to ensure data quality prior to our analysis. We followed the approach of Weitzel and Codding64 in removing all dates which are (1) duplicated in the database (based on the laboratory number), (2) marked as “anomalous” by the uploader, and (3) reported to originate from paleobiological or geological (i.e., nonanthropogenic) contexts. Additionally, we removed radiocarbon dates (n = 54) derived from megafauna skeletal material that originated from sites where convincing kill/scavenging associations could not be made32 and thus would not relate to past human activity. This left us with a sample of 938 anthropogenic radiocarbon dates.

Although previous analyses using the CARD dataset for the Paleoindian period have taken the data at face value while recommending caution23, the contested nature of pre-Clovis (>~13.2 ka) occupations warrants further scrutiny. We thus conducted separate analyses using the CARD dates between 13.2 and 10.0 ka but for the pre-13.2 ka period, only dates from the few most widely accepted pre-Clovis sites that have produced substantial radiocarbon records; these include Paisley Caves, Meadowcroft Rockshelter and Page-Ladson65,66,67 (see Supplementary Note 1: Analysis of Vetted Pre-Clovis Dataset, Supplementary Tables 1 and 2, and Supplementary Figure 1). The results and conclusions based on these analyses do not differ substantively from those presented above using only the dates obtained from the CARD. As it has been argued that variation in economic complexity can influence per capita energy consumption and the production of dateable archeological materials (bone, charred wood, etc.), semi-independent of population size68, we observe that our archeological date samples were derived from human populations that shared broadly similar technologies, foraging adaptations, and levels of economic complexity.

As noted above, our radiocarbon datasets for extinct megafauna are based exclusively on directly dated specimens derived from paleontological contexts following the systematic evaluation of Grayson and Meltzer32. Since no such evaluation has been provided for the Bison radiocarbon record, our analysis of this taxon (Fig. 6) is necessarily based on the collective assemblage of dates derived from both archeological and paleontological contexts. Since the radiocarbon dates for Bison are, to a large degree, a subset of the broader human radiocarbon record, a bison−human correlation analysis would be inappropriate and, thus, was not conducted. The sample of direct dates for Bison (n = 89) was compiled from the CARD58, Shapiro et al.19, and McDonald60.

Constructing summed probability distributions

Using the radiocarbon date samples obtained for both megafauna and humans, we first calibrated all dates based on the IntCal13 calibration curve using the rcarbon package69 available for R. Calibrating dates using this package produces identical results to conventionally used software such as OxCal54,55,56,57. Using the rcarbon package, we then constructed SPDs for the contiguous US using a 200-year moving average from 10 to 20 ka for humans (n = 938) and all taxa represented by more than 20 dates for this period. These taxa include mammoth (n = 74), mastodon (n = 99), Shasta ground sloth (n = 40), horse (n = 25), and saber-toothed cat (n = 21). To examine megafaunal−human population relationships on more regional scales we derived separate SPDs for the Southwest and the Great Lakes regions as they are represented by the largest samples of megafaunal dates. We used a broadly encompassing definition of the Southwest here to include the states of California, Nevada, Utah, Arizona, and New Mexico; the Great Lakes region, for our purposes, includes all states and provinces that border any of the Great Lakes (i.e., Ontario, New York, Pennsylvania, Ohio, Michigan, Indiana, Illinois, Wisconsin, and Minnesota). For these two regions, we separately examined SPDs for the megafauna taxa represented by more than 20 dates, namely, Shasta ground sloth (n = 29) and mammoth (n = 27) in the Southwest, and mastodon (n = 80) and mammoth (n = 23) in the Great Lakes. For the period of greatest interest here (15.0−11.0 ka), date density (dates/year) ranges from 0.158 to 0.002 for the different regional megafauna and human date records with a mean density value of 0.023. These values fall within the range of most published analyses using SPDs to reconstruct population trends23,24,25,26,27,28,54,55,56,57. Furthermore, Timpson et al.55 determined that the method used here to construct SPDs is robust to small sample sizes. Randomly subsetting a regional radiocarbon SPD from Europe, they found that the shape of the curve was generally maintained with sample sizes as low as 12 dates spanning 4000 years (13% of the 93 dates used to construct their original SPD). However, while the general patterns of significant deviation from the null model were maintained in this small subset of dates, the exact timing of these deviations is susceptible to minor changes due to the widened confidence intervals associated with smaller sample sizes. SPDs constructed from small samples of dates may still therefore record accurate patterns.

Significant SPD deviations from the null models

We followed the approach of Shennan et al.54,55,56,57 and compared each SPD to a null model to identify periods of time where the SPD deviates significantly from expected values. The null model to which we compared each SPD represents the trajectory each megafaunal taxon or human SPD followed prior to its decline. To create such a null model, we first used segmented regression analysis to fit paired quasi-Poisson family generalized linear models (GLMs) with a log link function to the SPD using the segmented package70. This analysis identified “breakpoints” at which the two GLMs were joined. These breakpoints represent the date at which the SPD shifts from an overall pattern of increasing to decreasing values. We next created a single quasi-Poisson family GLM with a log link that was based on the SPD values from 20.0 ka to the taxon’s specific breakpoint date, thereby representing only the increasing portion of the SPD. To generate the values of the null model itself, we allowed this GLM to predict values beyond the taxon’s breakpoint date based on the pre-breakpoint SPD trajectory. The resulting null model therefore reflects the pre-breakpoint (and pre-decline) exponential growth rates of each SPD, or the expected trend of the SPD should it never decline.

Then, using the modelTest function in the rcarbon package, we used a Monte-Carlo method to simulate 500 SPDs based on this null model and generate a confidence interval to which we compare the original SPD. This was done by taking the number of dates used to construct the original SPD and distributing them proportionately under the null model, using the shape of the null model to provide probability estimates for each year according to which dates were simulated. These simulated dates were then “uncalibrated” by passing them back through the IntCal13 curve: for each simulated calendar date, a random value was selected from within the calibration curve’s error range for that given year. Standard deviations associated with uncalibrated dates were next sampled from our original dataset of 14C dates with replacement and attached to each simulated date, ensuring that the simulated dates have comparable error estimates to the originals. The simulated 14C dates were then recalibrated using these sampled error values and used to construct an SPD. This SPD was created using a 200-year moving average applied to each year between 10.0 and 20.0 ka. Such a moving average helps to control for any small-scale artifacts of the calibration curve that would otherwise show through in the resulting SPD. The SPD was next locally transformed using Z-scores to standardize the probability distribution and make it comparable to others, which also has the effect of “de-trending” the resulting SPD so that any features of the original SPD and the calibration curve were removed. We repeated this process 500 times and a 95% confidence interval was generated around the exponential null model based on the Z-scores of each SPD, making it possible to identify time periods for which the observed SPD significantly deviates from the null model. Finally, a “global” P value expressing the overall deviation of the entire observed SPD from the null model was generated by calculating the proportion of simulated SPDs for which the total area lying outside the 95% confidence interval was as or more extreme than that of the observed SPD. These steps were then repeated for each taxon evaluated in each of the three regions, as well as for all taxa from all regions using taphonomically corrected dates following Surovell et al.71 (see Supplementary Note 2: Analysis of Taphonomically Corrected Data, Supplementary Table 3, and Supplementary Figure 2–5).

The utility of comparisons to an exponential null model

Shennan et al.54,55 emphasize that this approach, based on an exponential null model, provides two key benefits. First, using an exponential null model provides a hypothesis of unregulated growth against which to test an SPD: a null hypothesis in which the population growth rate never declines. Of course, real-world populations are rarely unregulated, making a logistic growth model a generally more realistic representation of population dynamics. However, applying a logistic growth model to these data present challenges. Logistic population growth models reach their top asymptote due to the effect of carrying capacity, which constrains a population from growing beyond a certain point typically dictated by food supply. It is not possible to derive the terminal Pleistocene carrying capacity for each of these megafauna taxa, making it problematic to assign such values in constructions of null models. We do not know the available food supply for each taxon, nor should we assume that each taxon had reached its carrying capacity prior to its declines towards extinction. Nonetheless, to ensure that our choice of exponential null models does not alter our results in any substantial way, we repeated all analyses using logistic null models (see Supplementary Note 3: Analysis Using Logistic Null Models and Supplementary Figures 6–8). These models were fit using maximum likelihood estimation of five-parameter logistic models based on Richard’s Equation. A maximum likelihood approach permits the estimation of the carrying capacity parameter (and all other parameters) based on the characteristics of the observed data, not any a priori expectation. Such values may not reflect real-world carrying capacity, for example, but are an empirically robust way of fitting such a model to these data in the absence of other information. Attempts to fit the more commonly known three-parameter logistic model typically failed as most of these SPDs behave exponentially prior to their breakpoint dates, not logistically. Many of the logistic models fit here are therefore either poor fits to these data or obtain marginally better fits at the cost of several additional degrees of freedom. For this reason, these null models frequently resemble their exponential counterparts because maximum likelihood estimation places the top asymptote (i.e. carrying capacity) far above the SPD in order to optimize the model fit to the available data.

Second, an exponential model increasing towards the present takes the same form as the taphonomic correction curve first proposed by Surovell et al.71 and revised by Williams25. Taphonomic loss of dateable organic materials increases with time according to an exponential function; therefore, deviations of an SPD from an exponential null model are, in effect, deviations from a taphonomic correction curve. Shennan et al.54,55,56,57 therefore suggest that dates for which SPD values are significantly greater than the exponential null model represent not only higher populations than would be expected given unregulated growth (i.e., “booms”), but also greater representation of dated carbon than would be expected according to a taphonomic correction curve. Conversely, when SPD values are significantly lower than the null model, this reflects lower populations than would be expected (i.e., “busts”) as well as a dearth of dated carbon relative to the amount expected given the taphonomic correction curve. However, because this approach simply provides an analog for the taphonomic correction curves of Surovell et al.71 and Williams25 against which to compare an SPD, and does not actually correct the SPD itself, we also repeat all of our analyses using taphonomically corrected dates, following the equation provided by Surovell et al.71. No substantive differences are apparent between the analyses we presented above and those based on the taphonomically corrected data (see Supplementary Note 2).

Potential sources of error

Bias in researcher submission of samples to date is also a potential source of error in both the CARD and megafaunal datasets and may be especially concerning in the context of the Paleoindian period in general (including Clovis, Folsom and related fluted point complexes) that is associated with heightened interest among archeologists. Yet we also anticipate equal or greater interest and bias in date submission for potential pre-Clovis sites. Collectively, we thus expect some oversampling for the Pleistocene-Holocene transition (~15.0–10.0 BP) period in general but see no mechanism to account for submission bias and oversampling to fluctuate in a patterned way within it. It is also unclear how researcher submission bias could account for the trends in the paleontological megafauna database; there would seem to be no less inherent interest in a mammoth skeleton (or a piece of Shasta ground sloth dung) from geological/stratigraphic contexts suggestive of the YD period, compared to those from contexts immediately preceding it. And as noted above, we excluded from our megafaunal dataset those dates derived from clear archeological contexts to avoid the inflation of megafauna numbers during the period they overlapped in time with people, while also only including one date per megafauna individual.

Paleoclimatic data

The paleoclimatic data presented in Figs. 2–4, 6 (and Supplementary Notes 1–3) are derived from the δ18O values from the NGRIP ice core72, the CO2 data from the EPICA Dome C, Antarctica, ice core46, and the insolation seasonality values at 60°N45.

Code availability

The R code used in this study is provided in Supplementary Data 2.