Climate and demography drive 7000 years of dietary change in the Central Andes

Abstract

Explaining the factors that influence past dietary variation is critically important for understanding changes in subsistence, health, and status in past societies; yet systematic studies comparing possible driving factors remain scarce. Here we compile the largest dataset of past diet derived from stable isotope δ¹³C‰ and δ¹⁵N‰ values in the Americas to quantitatively evaluate the impact of 7000 years of climatic and demographic change on dietary variation in the Central Andes. Specifically, we couple paleoclimatic data from a general circulation model with estimates of relative past population inferred from archaeologically derived radiocarbon dates to assess the influence of climate and population on spatiotemporal dietary variation using an ensemble machine learning model capable of accounting for interactions among predictors. Results reveal that climate and population strongly predict diet (80% of δ15N‰ and 66% of δ13C‰) and that Central Andean diets correlate much more strongly with local climatic conditions than regional population size, indicating that the past 7000 years of dietary change was influenced more by climatic than socio-demographic processes. Visually, the temporal pattern suggests decreasing dietary variation across elevation zones during the Late Horizon, raising the possibility that sociopolitical factors overrode the influence of local climatic conditions on diet during that time. The overall findings and approach establish a general framework for understanding the influence of local climate and demography on dietary change across human history.

Figure 1

Study area and spatial climate variation. Left: (a) Map of locations from which individuals are derived. Colors correspond to the elevation zone. Right: Averaged climate value for each grid cell from the past 7000 years for the study area. Climate variables are (b) temperature (°C), (c) temperature seasonality (standard deviation in °C * 100), (d) precipitation (mm/day), and e) precipitation seasonality (standard deviation in mm/day * 100) generated from TraCE21ka simulations^, through the PaleoView Paleoenvironmental Reconstruction Tool. Maps were produced in the R statistical environment using the mapdata and raster packages with shapefiles obtained from the NOAA Global self-consistent, hierarchical, high-resolution geography database.

Figure 2

Time series of climate and demographic variables. (a–d) z-scores of temporal climate deviation from the mean for the past 7000 years. Each grey line represents the z-score variation of an individual grid cell (see Fig. 1). Colored lines represent the moving average z-score for all raster grid cells combined. In order, the variables are (a) mean temperature (°C), (b) temperature seasonality (sd °C * 100), (c) mean precipitation (mm/day), and (d) precipitation seasonality (sd mm/day * 100). (e–j) Coastal, mid-elevation, and highland KDE estimates documenting relative variation in population size over time for each elevation category employing uncorrected (left) and taphonomically corrected (right) KDEs. For analytical purposes, individuals receive climate value estimates from the time series of the grid cell at their spatial coordinates, not the moving average central tendency line.

Figure 3

Generalized additive model time series. GAM regression results for both (a) δ¹⁵N‰ and (b) δ¹³C‰ across the three elevation categories. Dots represent the observed individuals plotted at their median date determined via resampling of their date range 10,000 times (weighted by calibrated radiocarbon probability for directly dated individuals). Solid lines are the central tendencies for each elevation category with shaded 95% confidence intervals. Lines were generated by creating 10,000 GAMs per elevation zone and isotope (60,000 total) where, for each GAM, each individual received a single year as their date, sampled via weighted sampling of their date range. The central tendency line is the mean fit line of the 10,000 GAMs for each sample, with 95% CIs created from the mean standard error and standard deviation of the 10,000 GAMs for each sample. Time periods follow the chronological periods from Moseley: Preceramic (PC), Initial Period (IP), Early Horizon (EH), Early Intermediate Period (EIP), Middle Horizon (MH), Late Intermediate Period (LIP), Late Horizon (LH). Note: This is a time series only and does not incorporate space. Some of the variation within elevation zones relates to spatial variation which is addressed in the RF analyses.

Figure 4

Cumulative climate and demographic effect sizes. Cumulative climate and demographic effect sizes (percent of the per mil amount of variation explained) within each elevation category for δ¹⁵N‰ (left column) and δ¹³C‰ (right column). (a, b) The cumulative effect sizes employing the uncorrected KDEs. (c, d) The cumulative effect sizes employing the taphonomically corrected KDEs. Regardless of KDE used for demographic estimates, climate change consistently has a larger effect than demography.