This is the submitted version of the journal article: Zhu, Dan; Galbraith, Eric; Reyes-García, Victoria; [et al.]. «Global hunter- gatherer population densities constrained by influence of seasonality on diet composition». Nature ecology & evolution, Vol. 5, Issue 11 (November 2021), p. 1536-1545. DOI 10.1038/s41559-021-01548-3 This version is available at https://ddd.uab.cat/record/276627 under the terms of the license 1 Global hunter-gatherer population densities constrained by influence of 2 seasonality on diet composition 3 Authors: Dan Zhu1,2*, Eric D. Galbraith3,2,4, Victoria Reyes-García2,4, Philippe Ciais5,6 4 Affiliations: 1 5 Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, 6 Peking University, Beijing, China 2 7 Institut de Ciència i Tecnologia Ambientals (ICTA-UAB), Universitat Autònoma de Barcelona, 8 08193 Barcelona, Spain 3 9 Department of Earth and Planetary Sciences, McGill University, Montreal, QC H3A 0G4, 10 Canada 4 11 ICREA (Catalan Institution for Research and Advanced Studies), 08010 Barcelona, Spain 5 12 Laboratoire des Sciences du Climat et de l’Environnement, IPSL-LSCE, CEA-CNRS-UVSQ- 13 UPSACLAY, 91191 Gif sur Yvette, France 6 14 Climate and Atmosphere Research Center (CARE-C) The Cyprus Institute 20 Konstantinou 15 Kavafi Street, 2121, Nicosia, Cyprus 16 17 *Correspondence to: zhudan@pku.edu.cn 18 19 1 20 Abstract: 21 22 The dependence of hunter-gatherers on local net primary production (NPP) to provide food 23 played a major role in shaping long-term human population dynamics. Observations of 24 contemporary hunter-gatherers have shown an overall correlation between population density 25 and annual NPP, but with a thousand-fold variation in population density per unit NPP that 26 remains unexplained. Here we build a process-based hunter-gatherer population model 27 embedded within a global terrestrial biosphere model, which explicitly addresses the extraction 28 of NPP through dynamically-allocated hunting and gathering activities. The emergent results 29 reveal a strong, previously unrecognized effect of seasonality on population density via diet 30 composition, whereby hunter-gatherers consume high fractions of meat in regions where 31 growing seasons are short, leading to greatly reduced population density due to trophic 32 inefficiency. This seasonal carnivory bottleneck largely explains the wide variation in population 33 density per unit NPP, and questions the prevailing usage of annual NPP as the proxy of carrying 34 capacity for ancient humans. Our process-based approach has the potential to greatly refine our 35 understanding of dynamical responses of ancient human populations to past environmental 36 changes. 37 38 39 2 40 Introduction 41 Hunter-gatherer populations, that subsist on hunting, gathering and fishing, rely intimately on the 42 biotic fabrics of their local environments1–3. Although only a handful of hunter-gatherer 43 populations remain, over 300 societies of contemporary hunter-gatherers have been documented 44 sporadically over the past two centuries1. Due to the assumed similarities between their lifestyles 45 and those of our foraging ancestors, contemporary hunter-gatherers have provided many insights 46 on our species’ past4–8, informing reconstructions of population changes on evolutionary 47 timescales5,6, the timing and rates of global human dispersals4,7, and the drivers behind the origin 48 of agriculture8. Yet, these contemporary societies do not provide direct analogues of ancient 49 foragers, as they had more complex technologies9 and had experienced interactions with and 50 pressures from neighbouring agricultural and industrial societies10,11, such as acquiring 51 supplemental agricultural food2,12 or encountering novel pathogens near the time of 52 documentation, which may have modulated the population density-ecosystem relationships. 53 Thus, many of the observational data may provide distorted views on the pre-existing state, 54 casting doubt on statistical models and demographic parameters that are directly fitted to the 55 contemporary data5–8. Taking full advantage of the insights from contemporary hunter-gatherers 56 requires a mechanistic, process-oriented understanding of how environmental factors influence 57 the distribution and abundance of hunter-gatherer populations. 58 Because hunter-gatherers acquire food directly from their surrounding environment, it has long 59 been thought that their population density should be closely linked to the productivity of their 60 local ecosystems1,2, and numerous studies have used net primary production (NPP) as the main 61 predictor for ancient hunter-gatherer density changes4,7,13. However, although contemporary 62 hunter-gatherer population density is positively correlated with NPP, NPP alone explains less 63 than 30% of the variability of population density (Fig. 1a). The fraction of NPP consumed by 64 hunter-gatherers (ΦNPP, see Methods), which also indicates population density per unit NPP, 65 varies by three orders of magnitude across the recorded groups (Fig. 1b and Supplementary 66 Figure 1e). A portion of this scatter could reflect the confounding historical influences, but it has 67 also been assumed that the edible proportion of primary production differs among biomes2. 68 However, our analysis does not show systematic differences in ΦNPP between biome categories, 69 except for a lower ΦNPP in boreal forest and a marginally significantly higher ΦNPP in 70 Mediterranean forest (Fig. 1b and Supplementary Figure 2b). Nor does ΦNPP show a consistent 71 change across the NPP gradient (Supplementary Figure 2c). The wide spread in ΦNPP thus 72 remains unexplained. The weak statistical dependence of population density on NPP has 73 prompted a recourse to other explanatory variables, including biodiversity3, pathogen stress3, 74 precipitation seasonality8, climate variability14, and social complexity-related variables such as 3 75 food storage-dependence14. It is however difficult to tease apart the causality among these inter- 76 correlated variables using statistical methods alone. 77 Here we build a mechanistic understanding of hunter-gatherer populations with a global, 78 activity-based process model for which population density is an emergent feature. Our model 79 operates within the framework of an Earth system model (ESM) (Fig. 2). During the past 80 decades, the land components of ESMs have integrated vegetation modules to resolve key 81 processes of biogeochemical cycling of carbon, energy and water15, and have recently begun to 82 incorporate large mammalian herbivores based on metabolic and demographic equations16,17. 83 Our model directly couples human population dynamics to such a land model17, simulating 84 human time allocation to hunting versus gathering in response to the interaction between humans 85 and food resources (see Methods). The model uses explicit formulations of daily carbon/energy 86 flows among vegetation, herbivores and humans, the outcome of which determines human 87 reproduction and mortality rates. These fluxes depend on a time allocation scheme with two 88 simple assumptions: first, total foraging time increases or decreases depending on the level of 89 fullness of the previous day, and second, gathering vs. hunting time depends on the relative 90 abundance of plant vs. animal food and on an underlying a priori preference for meat. The model 91 is resolved on a daily time-step, thus capturing the seasonal cycle in high detail. 92 93 Results 94 Modelled global distribution of hunter-gatherers 95 Simulated steady-state global population density and foraging times under present-day climate 96 are shown in Fig. 3. The model simulates high densities in most regions with high NPP and 97 closely matches the observed relationship of density versus NPP derived from contemporary 98 hunter-gatherer data (Fig. 3b). Simulated foraging times are generally 3.5-5.5 hours per person 99 per day throughout the world, close to the average of available observations from a few hunter- 100 gatherer groups2,18,19 (ca. 4 hr. day-1 ind.-1). In the model, slightly more time is spent hunting than 101 gathering in the tropics, whereas most of the foraging time is devoted to hunting at high latitudes 102 (Fig. 3cde). 103 The simulated population density agrees reasonably well with the ethnographic data (Pearson 104 correlation coefficient r=0.58, Supplementary Figure 3ab), and there is also a strong correlation 105 between simulated and ethnographically recorded diet compositions, though less of the observed 106 variance is explained (r=0.40 for the meat fraction of the diet, Supplementary Figure 3cd). 107 Model-data discrepancies may be partly due to the lack of fishing in the model, which is a major 108 source of food in many coastal societies (Supplementary Figure 1d). In addition, the model treats 4 109 all plant or animal food equally, thus ignoring the wide variety such as the different energy 110 content, protein richness, and relative abundance of each food item within both broad 111 categories2,20–22. Any bias in the vegetation and herbivore dynamics model17 could also 112 propagate to the human model, contributing to the model-data mismatch (note that the S2 113 experiment has corrected part of the biases in the vegetation inputs to the FORGE model, see 114 Methods). Furthermore, discrepancies could arise from climate shifts and non-steady state 115 factors given that the model simulates equilibrium densities forced by present-day climate (while 116 the ethnographic records were collected during the past two hundred years1) whereas some of the 117 populations were under positive growth rate when being studied23. There are also significant 118 uncertainties in the diet composition data, as the estimation methods for the amounts of gathered, 119 hunted, and fished food are often inconsistent, sometimes measured in different units (weight or 120 calorie) and occasionally based on the ethnographers’ impressions2, which may further add to the 121 model-data scatter. 122 The modelled global total hunter-gatherer population is 17 million, which is at the high end of 123 the estimates of prehistoric (pre-agriculture) population derived by extrapolating national 124 historical records, which range from 1-20 million24. This can be partly explained by the fact that 125 our model was calibrated with data from contemporary hunter-gatherers, who presumably had 126 more advanced technologies and access to non-local foods than early-Holocene foragers25 (see 127 Supplementary Discussion 2 for the test runs regarding the technology-dependent foraging 128 efficiency parameter). Geographical differences also contribute to our estimate being higher than 129 previous upscalings of contemporary hunter-gatherer populations, which include regression 130 models (10~12 million, refs.3,26) and a proportional projection method (7 million, ref.1). For the 131 former, the difference is primarily attributable to a higher density in Africa and southern Asia in 132 our model (Supplementary Figure 4), where the populations were suppressed in ref.3 due to a 133 high pathogen stress as a predictor in their regression (see Supplementary Discussion 1 for more 134 discussion about pathogen stress). For the latter, we find that the global estimate with the same 135 projection method but a different upscaling biome map ranges from 10 to 17 million, contingent 136 on how an average population density is calculated for each biome (see Supplementary 137 Discussion 3). Therefore, we consider our result of 17 million, simulated by a process-based 138 model, to be a reasonable, independent estimate of potential global hunter-gatherer population 139 size under modern climate with advanced foraging technology, assuming foraging as the only 140 subsistence type. 141 5 142 Linkage among seasonality, diet composition and population density 143 As in the observations, simulated ΦNPP (fraction of NPP consumed by hunter-gatherers) displays 144 a huge spatial variation with a thousand-fold range (Fig. 4b). The spatial distribution of ΦNPP 145 follows an inverse relationship with the meat fraction of the diet (compare Figs. 4a and 4b). The 146 negative correlation is clear in the scatter plot, both for grid cells where the observational data 147 are located (Fig. 4d) and for all populated grid cells in the model (Supplementary Figure 5), 148 showing a drop of ΦNPP by almost an order of magnitude as meat fraction increases by 60%. 149 Remarkably, this emergent property of the model is corroborated by ethnographic data, which 150 falls over a comparable range (compare Figs. 4c and 4d). Indeed, in both the observations and 151 the model, a higher meat fraction of the diet is associated with a lower population density under 152 the same level of NPP (Supplementary Figure 6). It is possible that the use of aquatic resources, 153 which are not included explicitly in the model, could increase the density of some populations 154 relative to terrestrial NPP (and thus ΦNPP) and thus confound the model-observation comparison. 155 However, we checked that the observed correlation between ΦNPP and meat fraction holds when 156 the groups with fish-rich diets are excluded (Supplementary Figure 7). Therefore, although 157 aquatic resources are important for hunter-gatherers in some locations (Supplementary Figure 158 1d), the association between ΦNPP and meat fraction appears to be a robust feature of the 159 observations. Notably, even the groups with intermediate (30%≤fishing<50%) or high (≥50%) 160 fishing percentages display significant negative correlations between ΦNPP and meat fraction as 161 well (Supplementary Figure 7), lending confidence to this overall trend. 162 The meat fraction of the modelled hunter-gatherer diet is mainly driven by a scarcity of 163 vegetation edible by humans during winter or dry seasons in the model. In regions with relatively 164 small seasonal variations in NPP, foragers can subsist mainly on gathering throughout the year 165 due to a substantially higher abundance of plant food than animal food at equilibrium (Fig. 5a). 166 By contrast, in regions with long non-growing seasons, edible plant biomass becomes depleted 167 and foragers depend on hunting to subsist through the season of scarcity. As a result, the 168 production rate of animals imposes a strong constraint on human population density in regions 169 with short growing seasons, regardless of how abundant plant food may be during the growing 170 season (Fig. 5b). Since relatively little energy flows from primary to secondary production due to 171 trophic inefficiency (Fig. 5), the ΦNPP is expected to be lower for humans subsisting more on 172 animals. Therefore, the apparent negative correlation between ΦNPP and meat fraction (Fig. 4d 173 and Supplementary Figure 5) is caused both directly by the effect of shifts in humans’ trophic 174 position, and indirectly by the covariation between meat fraction and seasonality, the latter being 175 able to influence ΦNPP through the limiting seasonal minima. 6 176 With the process-based model, we quantify the contribution of the two effects, by conducting 177 two experiments in which humans only subsist on either plants or animals (via allocating all 178 foraging time in either gathering or hunting, hereafter “onlyGather” and “onlyHunt”). The 179 difference in the resulted ΦNPP between the two runs could be regarded as the direct impact of the 180 trophic position. In the results of onlyGather, humans only persist in regions with long growing 181 seasons (> ca.180 days), but the total population size in these populated areas is 9.5 times higher 182 than that in the onlyHunt experiment for the same area (Supplementary Figure 8). 183 Correspondingly, ΦNPP is mostly higher in onlyGather than in onlyHunt (albeit the opposite 184 occurs in a few grid cells with relatively shorter growing seasons where a full reliance on hunting 185 slightly increases population density), and log10(ΦNPP) increases by 0.85 (or, ΦNPP increases by 186 seven times) on average in onlyGather compared to onlyHunt (Supplementary Figure 8e). 187 Therefore, the direct effect of diet composition via trophic energetics contributes about half to 188 the apparent slope of ΦNPP versus meat fraction. 189 The modelled contrast of seasonal cycles in diet composition between long versus short growing 190 season regions (Fig. 5) is qualitatively supported by observations in the Ache (in tropical 191 forest)20, and in the Hiwi (in savanna)27 and Bushmen groups (in Kalahari desert)28. The former, 192 living in places with relatively high NPP all year round, exhibits little seasonal variance in both 193 meat and plant intake, whereas the latter two, living in highly seasonal environments, exhibit a 194 highly varied share of plant food across seasons. An in-depth quantitative evaluation is however 195 difficult due to the limited number of hunter-gatherer groups with detailed, seasonally based diet 196 records, and because the model cannot capture the highly diverse phenology of various plant and 197 animal species at the local ecosystems. 198 Given the strong role of meat fraction in controlling population density predicted by the process- 199 based model, we carried out a statistical analysis of contemporary hunter-gatherer data in order 200 to test if meat fraction emerges as a significant predictor (see Supplementary Discussion 1). Our 201 analysis shows that, indeed, a multivariable stepwise linear regression with nine empirical 202 predictors (NPP, tree and grass coverage, biodiversity, growing season length, percentages of 203 food derived from hunting, gathering and fishing, and absolute latitude) identifies the hunting 204 percentage as the strongest explanatory variable for population density, followed by NPP and 205 grass cover, while the other variables are non-significant. Hunting percentage itself can be 206 explained by growing season length, followed by grass coverage (see Supplementary Discussion 207 1). Whereas it has been previously noted that the fraction of animal food tends to increase at 208 higher latitudes25,29, our result suggests that it is mainly the short growing season that explains 209 this gradient, as the effect of latitude itself is non-significant when growing season length is 210 considered. A structural equation model, which considers both direct and indirect effects of the 7 211 predictors, also confirms the dominant effect of hunting percentage on population density in the 212 observational dataset, and supports the link that seasonality influences population density via the 213 hunting percentage (Supplementary Figure 9 and Supplementary Table 1). Relationships 214 between ΦNPP and other variables, apart from hunting percentage (Fig. 4c), were also examined 215 (Supplementary Figure 10), showing a significant positive correlation between ΦNPP and growing 216 season length (r=0.37, p<0.001), which however becomes non-significant after controlling for 217 hunting percentage in the partial correlation analysis (Supplementary Table 2). 218 Discussion 219 This study built a process-based hunter-gatherer dynamics model, coupled to a global terrestrial 220 biosphere model, to explore the influence of environmental factors on hunter-gatherer population 221 density. The model explicitly simulates micro-scale processes including daily carbon/energy 222 flows among plants, animals and humans, combined with a dynamic time allocation for hunting 223 and gathering activities with simple assumptions. The emergent macroscopic relationships from 224 the model are well supported by ethnographic observations of contemporary hunter-gatherers, 225 revealing a causal mechanism whereby short growing seasons drive high fractions of meat in the 226 diet, leading to greatly reduced population density per unit NPP largely as a result of trophic 227 inefficiency. 228 The empirical negative association between hunting fraction and population density has been 229 reported or implied in previous studies14,30, but it was assumed that population density was the 230 causal factor, rather than vice versa: a higher population density was suggested to lead to less 231 frequent residential moves, more food storage, and thus greater dependence on the stationary and 232 higher productive food of plants and coastal resources14. Our process-based coupled human- 233 ecosystem model suggests a more parsimonious explanation, given that it reproduces the 234 observed relationships as a straightforward consequence of bioenergetics. Reinforcing 235 mechanisms between demography and diet via cultural factors may still be relevant, as suggested 236 by ref.14 and by anthropological studies on diet from both ecological/biological and sociocultural 237 perspectives31. In particular, food storage and sharing have received much attention, due to their 238 close association with societal characteristics and ability to mitigate food scarcity at daily to 239 seasonal timescales32–34. It may be interesting to include food storage and sharing in our model 240 framework in the future to quantitatively explore the interplay between environment, 241 demography, diet, and sociocultural aspects. 242 The integral links between seasonality, diet composition and population density have strong 243 implications for studies that explore ancient human populations. Although seasonal cycles in the 244 availability of different types of food have left traces in fossil teeth of extinct hominins35 and 8 245 been detected in the gut microbiome composition of the Hadza foragers in East Africa36, 246 seasonality has not been included in prior modelling studies4,6,7,13 that used only annual NPP 247 and/or mean paleoclimate variables, which would have biased expectations for the spatio- 248 temporal pattern of early human distribution and migration. Furthermore, the lengthening of 249 growing seasons following the last deglaciation, combined with the end-Pleistocene megafaunal 250 extinctions that occurred in different regions at different timing and rates37, might have driven 251 long-term changes in diet composition towards less meat-dependence. Our results suggest that a 252 shift from meat-dominated to plant-dominated diet could have boosted population much more 253 than implied by NPP changes alone. Given that demography might be at the core of cultural 254 evolution38, which further shaped human societies8,39, it is crucial to study the trajectory of past 255 population growth and its drivers, under the changing paleoclimate, vegetation and animal 256 distributions. 257 The model developed here represents a new type of computational modeling for ancient human- 258 ecology studies, breaking from a tradition40 that has mainly included niche (species distribution) 259 models based on statistical methods and without human feedbacks on the environment, and 260 agent-based models which resolve individual behaviors and interactions but require very detailed 261 local information that is challenging to assess and have only been applied at a local/community 262 scale41,42. By incorporating hunter-gatherers in a realistic, interactive and dynamical global 263 environment, the ESM framework helps to clarify how short growing seasons force humans to 264 rely on meat-rich diets, reducing ΦNPP and thus the population density per unit NPP. Additional 265 factors beyond those resolved in our current model, such as technological transitions25,43 and 266 long-range migrations, are sure to have played important roles in determining changes in hunter- 267 gatherer population abundance and diet through time. With further development, the approach 268 has the potential to provide more highly resolved pictures of the complex, multifaceted 269 interactions between our ancestors and the Earth system. 9 270 Methods 271 272 Contemporary hunter-gatherer data 273 The hunter-gatherer population density and coordinates were acquired from ref.3, which 274 combined the ethnographic datasets from refs.1,2,44 (357 data points in total). For the diet 275 composition (percentages of hunting, gathering, and fishing), we combined the two datasets of 276 subsistence from Kelly2 and Binford1,45 (340 data points in total), which indicate proportions of 277 food derived from terrestrial animals, terrestrial plants, and aquatic resources1, in terms of weight 278 or calorie2. Average of the two datasets were used if a population is present in both sources. 279 To explore the potential impact of environmental variables on population density and diets, 280 we extracted the values of the following variables from their global maps: annual NPP from the 281 satellite-derived MOD17A3 product46 (averaged over 2001-2010); fractional coverage of trees 282 and grasses from a satellite-derived vegetation cover product47; biodiversity index from ref.3, 283 which was the average of scaled richness of animals (mammals + birds) and vascular plants; 284 growing season length (GSL) calculated as the days in a year when daily gross primary 285 production (GPP) exceeds 20% of the maximum daily GPP for each grid cell, wherein the daily 286 GPP was from an upscaled global map based on FLUXNET tower sites48 (different thresholds 287 were tested but made little difference, e.g. GSL defined by a threshold of 10% and 20% are 288 highly correlated, r=0.93). Due to slight differences of the land-sea mask of these global land- 289 only products, a few coastal populations were located in grid cells with a non-valid value; in 290 such case, the value of the nearest land pixel was used. These potential environmental predictors 291 are used in statistical analyses of the ethnographic data as detailed in Supplementary Discussion 292 1. 293 Unlike Tallavaara et al.3, we did not use the climate-based Miami model to calculate NPP 294 but instead used MODIS NPP, because the Miami model was shown to overestimate NPP in the 295 tropics since it only considers annual mean temperature and precipitation and misses important 296 factors like nutrient limitation in tropical forests49. It was argued in ref.3 that the satellite-based 297 NPP was unfavourable because of the recent human appropriation of NPP. However, we 298 consider it a minor problem here because the majority of the sites have very small coverages of 299 cropland, according to a satellite product50 (242 out of the 357 sites have <5% cropland cover, 300 305 sites <20%). 301 We calculated the fraction of NPP consumed by hunter-gatherers (ΦNPP), which represents 302 the proportion of total energy available to heterotrophs that flows into hunter-gatherers, as 303 follows: per capita daily energy consumption, averaged over eight different groups2, is 2480 kcal 304 ind.-1 day-1, which equals to 174 kgC ind.-1 yr-1, using the conversion coefficients of 9.8 MJ 305 kgDM-1 (see the text below Eq. 4) and 0.45 gC gDM-1. This same value is multiplied by 306 population density and divided by NPP to derive ΦNPP for each site. 307 A potential systematic bias in the ethnographic diet composition data is worth attention. In 308 this study we regard the recorded hunting fractions as representing meat fractions of the diet. 309 However, the hunting fractions are likely to underestimate the true consumption of animals due 310 to the fact that early observers had devalued women’s contribution in meat provision, including 311 small animals, eggs and insects. Such contributions had been often categorized as “gathering”, or 312 simply ignored, as the focus of “hunting” had been on relatively large, mobile prey. This bias 313 could be large for groups in arid regions, e.g., the Western Desert in Australia51,52, for which the 10 314 recorded hunting percentage is about 30% while the true meat fractions are up to 80% when 315 women’s contribution to hunting is included51. The bias appears to be lower for groups in the 316 tropics like Amazon Basin and Congo Basin53, where recent field observations show less than 317 10% of total hunted food from women53. Unfortunately, revisiting and verifying the secondary 318 data on diet composition is difficult due to the disappearance of the foraging lifestyles for most 319 of the recorded hunter-gatherer groups. We therefore conducted a sensitivity test to see how the 320 observed negative relationship between ΦNPP and hunting fraction (Fig. 4c) might be influenced 321 by expected uncertainties and biases in the ethnographic dietary records. As detailed in 322 Supplementary Discussion 1, we found that the relationship between ΦNPP and hunting fraction is 323 robust to the expected biases, with only a minor impact on the slope of the correlation over the 324 most probable range of bias. 325 326 Description of the process-based hunter-gatherer dynamics model 327 Model overview 328 We designed the hunter-gatherer dynamics model, FORGE (FORager dynamics in Global 329 Ecosystem model), in the framework of Earth system models. It is coupled to the ORCHIDEE 330 (Organising Carbon and Hydrology In Dynamic Ecosystems) terrestrial biosphere model54, 331 which is the land component of the IPSL (Institut Pierre Simon Laplace) ESM. ORCHIDEE 332 simulates the energy and water balance, vegetation dynamics, and carbon cycle of land 333 ecosystems. Inputs for ORCHIDEE include meteorological variables (air temperature, 334 precipitation, incoming short and long wave radiation, wind speed, air humidity, and air 335 pressure), atmospheric CO2 concentration, land cover, topography and soil texture maps. Spatial 336 resolution and domain of each simulation are customized depending on the inputted climate 337 forcing. In the model, each grid cell is occupied by a suite of plant functional types (PFTs), with 338 their fractional covers adding up to one. ORCHIDEE has recently incorporated a module of large 339 mammalian herbivores that simulates animal population density driven by vegetation and 340 climate17. The main revision of the herbivore module compared to ref.17 is the inclusion of 341 browsers (herbivores that eat woody plants) whereas the previous version only calculated 342 grazers. The same parameterizations are used for browsers as grazers except that browsers feed 343 on the leaf and fruit compartments of tree PFTs, assuming 10% of these biomass is reachable to 344 browsers. The summation of grazer and browser density in each grid cell provides the potential 345 animal food density for humans. 346 The new FORGE model (Fig. 2) couples hunter-gatherers with vegetation and herbivores in 347 each grid cell via daily foraging activities and resultant carbon/energy flows, which then updates 348 an energetic reserve (body fat) that determines birth rate and mortality rate of human 349 populations. The human foraging activities are formulated based on a time allocation scheme. 350 The feedback of human hunting on herbivore dynamics is taken into account, while the impact of 351 humans on vegetation growth is neglected in the current model, but the model infrastructure 352 leaves room for future developments regarding human-induced environmental changes (e.g. 353 changes in land cover, fire regimes, etc.) and subsequent impacts on the ecosystems. Major 354 simplifications of FORGE include that i) only an average human (body weight = 50 kg ind-1) is 355 represented; ii) fishing is not implemented in the current land-only model, leading to a potential 356 underestimation in population density of the groups for which marine resources provide much 357 food, mostly in coastal regions (Supplementary Figure 1); iii) human migration55 across grid 11 358 cells is not accounted for, although the fact that the ethnographically recorded average migration 359 distance per move and the total distance moved between camps in a year are ca. 27 km and 280 360 km, respectively1, means that most hunter-gatherers would not move over more than the 361 dimensions of one grid cell (on the order of 104 km2 each) in an annual cycle; and (iv) the lack of 362 representation of small animals in the model, including birds, reptiles, insects, and small 363 mammals (<5~10 kg ind.-1), could lead to underestimations in available animal food and meat 364 fraction of the diet, which coincides with the potential systematic underestimations in the 365 ethnographic hunting fraction data mentioned above. Detailed formulations of FORGE are 366 described below. 367 Daily intake 368 The model considers two types of subsistence activity, hunting and gathering. The 369 contribution of each activity to the diet depends on the abundance of animal biomass and edible 370 vegetation biomass, the time spent in hunting and gathering, as well as the technology-dependent 371 efficiency of both activities: 372 = 373 = (1) 374 where IV and IA are daily dry matter intake of vegetation and meat (kgDM day-1 ind.-1); BV 375 and BA are biomass density of edible vegetation and animals (kgDM m-2); tg and th are the time 376 allocated to gathering and hunting (hr day-1); a is the searching area per person hour, fixed at 377 4000 m2 hr-1 ind.-1; eg and eh are the efficiencies of gathering and hunting, representing the 378 fraction of the corresponding biomass that is acquired when hunter-gatherers pass the searched 379 area. 380 In a linear form this equation would imply that the food resources and foragers are both 381 uniformly distributed within each grid cell and foragers encounter them randomly. In reality, 382 excursions are undertaken from localized camps, different food items are distributed in patches, 383 and foragers know how to target the higher energy-return items first2. To partly account for these 384 aspects of sub-grid heterogeneity, we assume eg and eh to decrease as Bv and Ba decrease relative 385 to human density, capturing the increasing difficulty of foraging as the best resources are 386 sequentially depleted. We model this with an asymptotic function: 387 = 388 = (2) 389 where BH is human biomass density (kg m-2), equal to PH (population density which is 390 prognostically simulated in the model, see Eq. 11) multiplied by mean body weight (WH, 50 kg 391 ind-1); emax is the maximum fraction of the edible plant/animal biomass that can be acquired 392 when hunter-gatherers pass the searched area, fixed at 1 and tested in a sensitivity analysis (see 393 below); c determines the patch depletion rate; the value of c is fixed at 100, and is tested in the 394 sensitivity analysis. Notably, the inclusion of Eq. 2 produces a dynamic similar to that found in 395 ecological studies of per capita prey consumption by wild predators, which often decreases with 396 predator density for a given prey density56. 12 397 Energy expenditure and energetic reserve 398 The per capita daily energy expenditure E (MJ day-1 ind.-1) is calculated as: 399 = + + (3) 400 where Ec is the energy expenditure excluding those spent during foraging, fixed as 8.37 MJ 401 day-1 ind.-1 (i.e. 2000 kcal day-1 ind.-1); Eg and Eh are the energy spent per hour gathering or 402 hunting, fixed as 1.25 MJ hr-1 ind.-1 (according to ref.57) for both parameters in the model. 403 Based on daily food intake and energy expenditure, per capita energetic reserve (mainly in 404 fat cells) is updated daily as below. This is the key variable that indicates health condition and is 405 used to calculate birth rate and mortality rate (see below). 406 = 407 = + (4) 408 where F (kg ind.-1) is per capita energetic reserve; m (MJ kg-1) is the conversion coefficient 409 between energy and fat mass, set to 39.3 if IE (anabolism) (ref.58). 410 The gross energy value of plants and animals is around 3800-4800 for different plants (ca. 5000 411 for seeds) and 4400-5600 for different animal taxa (ref.59, in the unit kcal kgDM-1). We thus 412 assume the same value for both food types for simplicity. Further considering the ca. 30% energy 413 loss in excreta, and that the metabolisable energy content is less than the gross energy value of 414 food, we finally set mV and mA to 9.8 MJ kgDM-1. It should be noted that the energy density of a 415 mammal’s carcass depends on its fat content22, thus a fixed mA value in the current model is a 416 simplification and does not account for variations across species and over time. 417 There is an upper limit of F (Fmax), assumed to be 30% of body weight, i.e. 15 kg ind.-1 in 418 the model. Accordingly, whenever F is about to exceed Fmax, daily intake will be reduced so that 419 F stays at Fmax. 420 Time allocation 421 In the current model, we only simulate the time spent in foraging, and leave all other 422 activities such as childcare, domestic maintenance and socializing to an aggregate ‘other’ 423 category. The hunting, gathering and total foraging time (tf=th+tg) are updated every day in the 424 following two steps. First, tf increases or decreases by 1 hr day-1 depending on the changes in the 425 energetic reserve F during the previous day: ( − 1) − 1, ∆ >0 ( ) > 0.5 426 ( )= ( − 1) + 1, 427 ∈[ , , , ] (5) 428 th and tg are then changed proportionately. Such reduction of foraging time when humans 429 have eaten plenty of food in the previous day could represent idle relaxation, a phenomenon 430 widely observed in many animals60 as well as hunter-gatherers61. An upper limit on foraging 431 time would also be driven by the need to do the other essential activities not modelled here. 432 Foraging time reported for dozens of hunter-gatherer groups, averaged between males and 433 females, range between 0.8-6.8 hr day-12,18,19. We therefore set tf,min=0.5 and tf,max=8 hr day-1. 13 434 Second, the allocation between hunting and gathering depends on the relative energy benefit 435 of the two activities, as well as a craving for meat, parameterized as an exponential decay 436 function of meat proportion in the diet of the previous day: 437 ℎ = 438 = ( ) 439 = 440 = 441 = 442 ℎ ∈ [0.05,0.95] (6) 443 where thfrac is the fraction of tf allocated to hunting. The additional meat-craving term, Cmeat, 444 is introduced to represent that hunter-gatherers generally have a preference for meat, which may 445 reflect nutritional needs (e.g. fat and essential amino acids)62 as well as cultural importance2. The 446 parameter q is fixed at 2.5 but also tested (see Supplementary Table 3). In cases when fmeat falls 447 below 10% (the minimum observed hunting+fishing fraction), thfrac is set to 0.95. th and tg are 448 then calculated as × ℎ and × (1 − ℎ ) respectively. 449 This time allocation scheme is highly simplified. It does not consider, for example, seasonal 450 variations of other necessary activities such as providing shelter, extreme weather that makes 451 foraging activities risky or unworthy (e.g. the Hiwi in tropical savannah don’t forage during the 452 middle of the day to avoid heat27), and additional social/cultural factors. Further improvements 453 and tests of the scheme are possible in the future, given more observational data on human time 454 allocation dynamics. Yet, it provides an unprecedented approach to include the fundamental limit 455 of time on the ability to capture existing biomass. 456 Updating food resources 457 Not all vegetation and herbivore biomass are edible or accessible to humans. We denote the 458 edible and/or accessible fraction of vegetation and herbivores as ΦV and ΦA, respectively, which 459 are assumed to be fixed parameters in the current model over the globe (ΦV=0.015, ΦV=0.1, see 460 sensitivity tests below). This is a simplification; for example, many fruits in tropical forests grow 461 high in the canopy, making them harder to reach than the case in short-vegetation biomes2. 462 Unfortunately, robust estimates of the edible/accessible fraction for major biomes worldwide are 463 unavailable, and therefore we do not vary ΦV and ΦA between grid cells in the current model, in 464 order to avoid arbitrarily introducing spatial relationships of human density versus NPP. Apart 465 from fruits, hunter-gatherers can also eat roots and tuber, which are from the belowground 466 primary production of some herbaceous plants. We thus also include part of the belowground 467 grass NPP as potential food sources. 468 The BV (biomass density of edible vegetation) is updated on a daily time-step: 469 = ( + )− − . . 470 = (7) 14 471 Where λ (day-1) is the daily turnover rate of the edible vegetation biomass; ffruit is the 472 fraction of NPP (including both tree and grass) allocated to fruits, set to 6.5% in the current 473 model according to a synthesis63 of fruit fall observations across tropical and temperate forests; 474 fbelow is the fraction of grass NPP allocated belowground, calculated as a function of mean annual 475 precipitation (MAP, mm yr-1) following the empirical equation from ref.64, which captures the 476 higher belowground allocations in dryer ecosystems. The lower limit of fbelow is set to 20%. The 477 inputted NPP and λ are either from ORCHIDEE outputs or from observation-derived values (see 478 “Model setup” section for the different runs). 479 The herbivore dynamics module is coupled once per year with the human dynamics module. 480 Detailed descriptions of the herbivore module can be found in ref.17. In summary, herbivore 481 population density (PA,tot, ind. km-2) is updated on a yearly time-step using a logistic equation 482 (Eq. 8 in ref.17). The herbivore birth and mortality rates depend on an animal fat reserve pool (kg 483 ind.-1) that is updated daily according to their energy intake through grazing or browsing. 484 Within the human module, BA is updated on a daily time-step: 485 =− (8) ∑ 486 At the end of each year, the accumulated number of hunted animals ( ) is subtracted × . 487 from the herbivore population density PA,tot. The decreased PA,tot is passed to the herbivore 488 module to calculate a new PA,tot for the next year. BA is then updated and used in the human 489 module for the new year: 490 = , × 0.28 (9) 491 where WA is the mean body weight of herbivores, set to 180 kg ind.-1 in the current model; 492 0.28 is the conversion between live weight and dry mass after excluding water and bones. 493 Note that the herbivore population dynamics is also affected by seasonality, with relatively 494 higher animal starvation mortalities in regions with longer winters/dry seasons17. However, 495 temporally, simulated meat density is relatively stable throughout the year in the model (see Fig. 496 5), because the animal population density is updated only once per year. The annual update is 497 justified by the fact that the animals represented have lifespans of ~25 years17, so that seasonal 498 fluctuations in population numbers are strongly damped. Unlike humans, herbivores can feed on 499 plant litter (dead grasses and fallen leaves/fruits in the model) to live through non-growing 500 seasons, providing vital food for hunter-gatherers to survive the long winters/dry seasons in the 501 model. 502 It should be emphasized that the modelled feedbacks of human foraging are very different 503 for vegetation and animal density. Gathering is assumed to have negligible impacts on primary 504 productivity, since the Bv that is harvested represents the edible part of plants, which generally 505 do not contribute to plant growth, and for which harvesting can even aid in reproduction (e.g. 506 through the spreading of seeds65). In contrast, hunting can directly reduce the annual production 507 of herbivores by reducing PA,tot. However this effect would be small at low rates of hunting, for 508 which competition among herbivores places the limit on PA,tot. 15 509 Human population dynamics 510 The population birth rate (rbirth, yr-1) is assumed to depend on the average body condition, 511 represented by the daily-varying energetic reserve F. Thus, rbirth is calculated every day and 512 averaged over the year to be used in Eq. (11), 513 = (10) ( ) 514 where ra = 0.1, rb = 15, rc = 0.5. This sigmoidal function is similar to the equation used to 515 calculate herbivore birth rate17, with the parameters modified so as to give a maximum birth rate 516 of 10%, close to the recorded highest crude birth rate during the past two centuries (ca. 6%, ref66) 517 (see sensitivity tests of the three parameters below). 518 For mortality, we consider two processes. First, a background mortality rate (Mb) which is 519 the inverse of lifespan, fixed as =0.0125 yr-1. Second, starvation-induced mortality (Ms, yr-1), 520 caused by the exhaustion of body energy storage as represented by the complete depletion of fat 521 reserves. The Ms is calculated using a similar method as in the herbivore module17. Namely, we 522 assume a normal distribution of body fat within the population, with a mean μ=F and a standard 523 deviation σ=0.125Fmax; the probability that fat mass falls below 0 is taken as the value of Ms. 524 Similar to rbirth, Ms is also calculated every day and averaged over the year. 525 Finally, the annual dynamics of the human population density are calculated as: 526 = −( + ) (11) 527 where PH (ind. km-2) is the human population density for each grid cell. The PH is initialized 528 as P0 =10-3 ind. km-2, a lower value than the minimum density ever recorded (0.002 ind. km-2, 529 ref.1,2). Whenever PH falls below P0, PH is reset to P0. 530 Once each year, the reproduction energy cost (313 kcal day-1, or 478 MJ for a year, 531 according to ref.67) is subtracted from the energetic reserve F in order to account for the energy 532 input to birth: 533 ∆ =− (12) 534 Model setup 535 To derive the input variables for FORGE, we firstly ran ORCHIDEE for the globe at 2° 536 spatial resolution and for 200 years as a spin-up to reach the equilibrium of vegetation 537 production and biomass under present-day climate conditions. The climate forcing was the CRU- 538 NCEP reanalysis dataset68, repeating the years 2001-2010. The atmospheric CO2 concentration 539 was fixed at 380 ppm (the average level during 2001-2010). For the land cover, although 540 ORCHIDEE can simulate vegetation distributions equilibrated under the given climate 541 conditions, we instead prescribed an observation-based PFT map in order to reduce the bias of 542 ORCHIDEE outputs that would propagate to the FORGE model. This PFT map was based on 543 ESA CCI Land Cover map v2.0.7b for 2010, which was further merged with the LUH2 dataset 544 to generate a pre-industrial PFT map (year 850) with minimal crop coverages. See detailed 545 descriptions of the PFT map at https://orchidas.lsce.ipsl.fr/dev/lccci/. 546 The last 10 years’ outputs from ORCHIDEE are used to provide the inputs to FORGE. 547 These include the inputs required by the herbivore module: fractional cover of PFTs (aggregated 16 548 into three types: grass, tree, bare ground), carbon influx rates (i.e. NPP allocated to the edible 549 plant tissues, and influx to the edible litter pool) and decay rates of the edible pools, and annual 550 mean temperature; as well as inputs required by the human module: NPP, decay rate of the fruit 551 compartment, and annual precipitation (used in Eq. 7). The input variables are at a daily time 552 step. FORGE was then run for 300 years to reach the equilibrium of population density of both 553 herbivores and hunter-gatherers. The last 50 years were averaged and presented as the results. 554 Supplementary Figure 11 shows the ORCHIDEE-simulated NPP, in comparison with 555 MODIS NPP. ORCHIDEE can reproduce the overall pattern of the satellite-derived NPP, with a 556 reasonable agreement with MODIS NPP at the sites (Pearson correlation coefficient r=0.67). 557 However, ORCHIDEE underestimated NPP in tropical forests and arid ecosystems including 558 some classical hunter-gatherer areas like the Great Basin in North America, Kalahari in southern 559 Africa, and interior Australia (Supplementary Figure 11). The underestimation of annual NPP in 560 the dry regions is partly due to an underestimation of growing season length by ORCHIDEE 561 (Supplementary Figure 12). 562 In order to test the impact of the bias of ORCHIDEE-simulated NPP, and to minimize this 563 bias that would propagate to the FORGE-simulated hunter-gatherer densities, we conducted three 564 sets of simulations, with different input files while identical FORGE model: 565 S0: inputs to FORGE are directly from ORCHIDEE outputs; 566 S1: compared to S0, the daily NPP is multiplied by a scaling factor so that the annual NPP 567 (average for 2001-2010) equals to MODIS NPP for each grid cell; 568 S2: compared to S0, the daily NPP is replaced by the observation-based values, that is, 569 annual MODIS NPP interpolated into daily values according to an observation-derived daily 570 GPP product48. Note that the herbivore module in FORGE separates grazers feeding on grass 571 NPP and browsers feeding on tree NPP, and here we assume the same NPP values (per unit PFT 572 area) for tree and grass PFTs in the same grid cell. Besides, the ORCHIDEE-simulated daily 573 decay rates of edible plant tissues (λ) are replaced by a constant of 0.01 day-1 during growing 574 season (when NPPdaily>NPPdaily,max) and 0.05 day-1 during non-growing season (when 575 NPPdaily<=NPPdaily,max), in order to avoid an inconsistent timing of senescence (a period of high 576 decay rates of leaves and fruits) with the observed seasonality in NPP. The λ values are chosen 577 so that the integrated annual decay rates are generally comparable to those simulated by 578 ORCHIDEE (used in S0 and S1). 579 580 Modelled population densities and their relationships between NPP and meat fraction of the 581 diet from the S0 and S1 runs are shown in Supplementary Figures 13 and 14. Among S0, S1, and 582 S2, the emergent relationships between population density and NPP, and ΦNPP and meat fraction, 583 are similar (Supplementary Figure 14), except for a steeper slope of log-population density 584 versus NPP in S0 due to the lack of data points under the high end of NPP (as ORCHIDEE 585 underestimates NPP in tropical forests, see Supplementary Figure 11). Regarding the global 586 patterns of population density (Supplementary Figure 13), however, both S0 and S1 predict an 587 unrealistic absence of hunter-gatherers in many arid ecosystems including those in western 588 Plains in North America, Kalahari, and interior Australia, the classical hunter-gatherer territories. 589 These regions are populated in the S2 results (Fig. 3a). The difference between S1 and S2 (for 590 which the annual total NPP in the inputs are identical) therefore indicates the critical role of a 591 substantially long growing season length in sustaining hunter-gatherers in the less productive 592 ecosystems. Considering a prevailing underestimation in growing season length for the arid 17 593 ecosystems in state-of-the-art terrestrial biosphere models (Supplementary Figure 15), our results 594 highlight the need for a more realistic simulation of plant phenology in order to improve the 595 simulation of contemporary and ancient hunter-gatherers. Unless otherwise specified, the model 596 results shown are from the S2 experiment. 597 598 Sensitivity tests of parameters 599 600 We consider several parameters in FORGE, listed in Supplementary Table 3, to be highly 601 uncertain. We therefore conducted sensitivity tests using Sobol’s method69 (a variance-based 602 sensitivity analysis), which decomposes the variance of model output into fractions that can be 603 attributed to different inputs. Sobol’s method has the advantage of measuring sensitivity across 604 the whole input space, as well as accounting for nonlinearity and parameter interactions, so that 605 the total-order index (ST) indicates the importance of each parameter considering both its main 606 effect (first-order sensitivity index, S1) and the contribution of its interaction with other 607 parameters. 608 Sobol’s method requires thousands of runs, which is computationally expensive for global 609 experiments. We thus carried out the tests at two sites with contrasting characteristics. First, a 610 temperate forest with a long growing season and a low hunting fraction in the diet at 42N, 123W 611 in. Second, a boreal forest with a short growing season and a high hunting fraction at 56N, 69W 612 (Supplementary Figure 16). For each site, we conducted 9000 runs in which the parameters were 613 sampled within their ranges listed in Supplementary Table 3, using Saltelli’s sampling scheme 614 which is more efficient than random sampling. Supplementary Figure 17 shows how the varied 615 parameter values affect the corresponding equations. All the sensitivity tests were done using the 616 Python package “SALib” (https://salib.readthedocs.io/en/latest/index.html). 617 The resulting ST and S1 (Supplementary Figure 18) indicate different sensitivities to the 618 parameters between the two sites. For the hunting-dominated site, population density is most 619 sensitive to ΦA, followed by emax and c, whereas for the gathering-dominated site, population 620 density is sensitive to the four parameters ΦV, ΦA, emax, and c. The three parameters used to 621 calculate birth rate make negligible contributions to the variance in modelled population. Indeed, 622 birth rate impacts the time it takes to reach equilibrium population under a given set of 623 environmental conditions, but has a negligible effect on the equilibrium value itself. Nor is the 624 population density sensitive to the parameter q, which determines the meat craving response. 625 Supplementary Figure 19 further shows the quantitative response of model results to the 626 parameters. As expected, a higher ΦV (or ΦA) value increases population density at the 627 gathering-dominated (or hunting-dominated) site, while a higher c value decreases the population 628 at both sites. A higher emax value increases population density at both sites, but the beneficial 629 effect diminishes at high population densities for the hunting-dominated site, probably a result of 630 overhunting. Meat fraction in the diet is mainly determined by the relative values of ΦV and ΦA 631 (Supplementary Figure 18). Compared to population density and meat fraction, foraging time is 632 relatively more equally sensitive to all the tested parameters (Supplementary Figure 18), but its 633 results in the 9000 runs are relatively centralized across the tested parameter space for both sites 634 (Supplementary Figure 19). 635 For the standard configuration, we set ΦV=0.015, ΦA=0.1, emax=1, and c=200 so as to match 636 the average value of population density and meat fraction across the sites. Note that changing 18 637 these globally-constant parameter values increases or decreases population everywhere, but has 638 minor impacts on the modelled relationships with environmental variables (Supplementary 639 Figure 20). 19 640 Data availability 641 The contemporary hunter-gatherer data and environmental variables used in the analysis are 642 available in the supplementary materials. 643 Code availability 644 Source code (in Python) of the FORGE model and its output files (in NetCDF format) for this 645 study, including the three sets of global simulations (S0, S1, S2), are provided in a 646 Supplementary zipped file. The corresponding input files for the FORGE model are deposited in 647 the PANGAEA repository at https://doi.org/XXX. 648 649 Acknowledgments 650 D.Z. and E.D.G. acknowledge the financial support from the European Research Council under 651 the European Union’s Horizon 2020 Research and Innovation Programme (under Agreement 652 682602, to E.D.G.). D.Z. also acknowledges support from the National Natural Science 653 Foundation of China (grant number 41988101). V.R-G. acknowledges support from the 654 European Research Council under Agreement 771056. 655 Author contributions 656 D.Z. and E.D.G. conceived the study and model design. D.Z. built the model, performed the 657 analyses, and wrote the first draft. E.D.G. provided discussion and suggestions throughout the 658 process. V.R-G. and P.C. contributed to the interpretation of the results and writing of the 659 manuscript. 660 Competing interests 661 Authors declare no competing interests. 662 663 Figure legends 664 Fig. 1 | Contemporary hunter-gatherer population density vs. net primary production. (a) 665 Relationship between population density, from ethnographic records at 357 locations 666 (Supplementary Figure 1), and NPP, according to the MODIS satellite-derived product (see 667 Methods). The solid line gives the linear regression of log10(pop. density) vs. NPP ( = 9.6 × 668 10 − 1.53) with 95% confidence intervals shown in grey. (b) Violin plot for ΦNPP (fraction 669 of NPP consumed by hunter-gatherers, calculated as population density multiplied by a constant 670 intake rate and divided by NPP, see Methods) of populations located in different biomes 671 (Supplementary Figure 2). The white circles represent median values and the thick (or thin) 672 black bars the inter-quantile (or 5th–95th) ranges. Except for a lower ΦNPP in boreal forest and a 673 slightly higher ΦNPP in Mediterranean forest, ΦNPP does not differ significantly among the other 674 five biome types (ANOVA test p >0.05). 675 20 676 Fig. 2 | Schematic of FORGE (FORager dynamics in Global Ecosystem model), coupled to 677 the ORCHIDEE global terrestrial biosphere model. FORGE simulates daily energy (carbon) 678 intake from plants and animals and energy expenditure (black arrows); these energy in/out-fluxes 679 update an energetic reserve pool (body fat) every day, which indicates the average health state 680 and then impacts birth rate and mortality rate, determining changes in hunter-gatherer population 681 density every year. The daily intake rate depends on food abundance and on the time spent in 682 hunting and gathering, calculated with time allocation algorithms (see Methods for details). 683 684 Fig. 3 | Modelled population density and time allocation. (a) Map of hunter-gatherer 685 population density equilibrated under present-day climate (results in the S2 experiment, see 686 Methods). (b) Relationship between population density and NPP at the grid cells where the 687 observational data are located. The red solid line indicates linear regression of log10(pop. density) 688 vs. NPP ( = 11.7 × 10 − 1.78, N=334, r=0.72), with 95% confidence intervals shown in 689 red shading. For comparison, regression of the observations is shown in grey (same as those in 690 Fig. 1a). (c-e) Modelled time spent in foraging (c), separated into gathering (d) and hunting (e). 691 The grid cells where modelled population density is less than 0.2 ind./100km2 are shown in white 692 in the maps, and excluded in (b). 693 694 Fig. 4 | Meat fraction of the diet controls ΦNPP. (a,b) Modelled fraction of hunted food in the 695 diet (a) and ΦNPP (fraction of NPP consumed by hunter-gatherers) (b). (c,d) Relationship 696 between ΦNPP and meat fraction of the diet according to ethnographic records (c) and model 697 results (d). The black and red solid lines give the linear regressions of log10(ΦNPP) vs. meat 698 fraction from observations (c) ( = −1.93 − 3.68) and the model (c) ( = −1.78 − 3.68, 699 N=334, r= -0.73), with 95% confidence intervals shown in shading. The grid cells where 700 modelled population density is less than 0.2 ind./100km2 are shown in white in the maps, and 701 excluded in (d). 702 703 Fig. 5 | Seasonality, diet composition and carbon flows. Model results were averaged over all 704 grid cells of either gathering-dominated (modelled annual mean hunted food <40%, see Fig. 4a) 705 or hunting-dominated (hunted food ≥40%) regions, at dynamic equilibrium. Mean population 706 density of the two regions is 0.19 and 0.04 ind./km2, and mean meat fraction of the diet is 22% 707 and 55%, respectively. To synchronize the seasonal cycles between the two hemispheres, outputs 708 of the southern hemisphere were shifted by 6 months. The numbers on the links indicate 709 percentages of annual NPP that flow to herbivores, humans, and litter (dead plants) and soil 21 710 pools. Note that for legibility, the link widths are not strictly proportional to the magnitude of 711 flows. Spring: MAM; Summer: JJA; Autumn: SON; Winter: DJF. 712 713 714 References 715 1. Binford, L. R. Constructing frames of reference : an analytical method for archaeological 716 theory building using hunter-gatherer and environmental data sets. (University of 717 California Press, 2001). 718 2. Kelly, R. L. The lifeways of hunter-gatherers : the foraging spectrum. (Cambridge 719 University Press, 2013). 720 3. Tallavaara, M., Eronen, J. T. & Luoto, M. Productivity, biodiversity, and pathogens 721 influence the global hunter-gatherer population density. Proc. Natl. Acad. Sci. 115, 1232– 722 1237 (2018). 723 4. Eriksson, A. et al. Late Pleistocene climate change and the global expansion of 724 anatomically modern humans. Proc. Natl. Acad. Sci. 109, 16089–16094 (2012). 725 5. Gurven, M. D. & Davison, R. J. Periodic catastrophes over human evolutionary history 726 are necessary to explain the forager population paradox. Proc. Natl. Acad. Sci. 201902406 727 (2019). doi:10.1073/pnas.1902406116 728 6. Tallavaara, M., Luoto, M., Korhonen, N., Järvinen, H. & Seppä, H. Human population 729 dynamics in Europe over the Last Glacial Maximum. Proc. Natl. Acad. Sci. 112, 8232– 730 8237 (2015). 731 7. Bradshaw, C. J. A. et al. Minimum founding populations for the first peopling of Sahul. 732 Nat. Ecol. Evol. (2019). doi:10.1038/s41559-019-0902-6 733 8. Kavanagh, P. H. et al. Hindcasting global population densities reveals forces enabling the 734 origin of agriculture. Nat. Hum. Behav. 2, 478–484 (2018). 735 9. Porter, C. C. & Marlowe, F. W. How marginal are forager habitats? J. Archaeol. Sci. 34, 736 59–68 (2007). 737 10. Reyes-García, V. & Pyhälä, A. Hunter-gatherers in a Changing World. (Springer 738 International Publishing, 2017). 739 11. Lee, R. B. & Daly, R. The Cambridge encyclopedia of hunters and gatherers. (Cambridge 740 University Press, 1999). 741 12. Kitanishi, K. Seasonal Changes in the Subsistence Activities and Food Intake of the Aka 742 Hunter-Gatherers in Northeastern Congo. Afr. Study Monogr. 16, 73–118 (1995). 743 13. Timmermann, A. & Friedrich, T. Late Pleistocene climate drivers of early human 744 migration. Nature 538, 92–95 (2016). 745 14. Keeley, L. H. Hunter-Gatherer Economic Complexity and “ Population Pressure ”: A 746 Cross-Cultural Analysis. 373411, (1988). 22 747 15. Fisher, J. B., Huntzinger, D. N., Schwalm, C. R. & Sitch, S. Modeling the Terrestrial 748 Biosphere. Annu. Rev. Environ. Resour. 39, 91–123 (2014). 749 16. Pachzelt, A., Forrest, M., Rammig, A., Higgins, S. I. & Hickler, T. Potential impact of 750 large ungulate grazers on African vegetation, carbon storage and fire regimes. Glob. Ecol. 751 Biogeogr. 24, 991–1002 (2015). 752 17. Zhu, D. et al. The large mean body size of mammalian herbivores explains the 753 productivity paradox during the Last Glacial Maximum. Nat. Ecol. Evol. 2, 640–649 754 (2018). 755 18. Dyble, M., Thorley, J., Page, A. E., Smith, D. & Migliano, A. B. Engagement in 756 agricultural work is associated with reduced leisure time among Agta hunter-gatherers. 757 Nat. Hum. Behav. (2019). doi:10.1038/s41562-019-0614-6 758 19. Hill, K., Kaplan, H., Hawkes, K. & Hurtado, A. M. Men’s time allocation to subsistence 759 work among the Ache of Eastern Paraguay. Hum. Ecol. 13, 29–47 (1985). 760 20. Hill, K., Hawkes, K., Hurtado, M. & Kaplan, H. Seasonal variance in the diet of Ache 761 hunter-gatherers in Eastern Paraguay. Hum. Ecol. 12, 101–135 (1984). 762 21. Marlowe, F. W. et al. Honey, Hadza, hunter-gatherers, and human evolution. J. Hum. 763 Evol. 71, 119–128 (2014). 764 22. Cordain, L. et al. Plant-animal subsistence ratios and macronutrient energy estimations in 765 worldwide hunter-gatherer diets. Am. J. Clin. Nutr. 71, 682–692 (2000). 766 23. Gurven, M. & Kaplan, H. Longevity Among Hunter- Gatherers: A Cross-Cultural 767 Examination. Popul. Dev. Rev. 33, 321–365 (2007). 768 24. Klein Goldewijk, K., Beusen, A. & Janssen, P. Long-term dynamic modeling of global 769 population and built-up area in a spatially explicit way: HYDE 3.1. The Holocene 20, 770 565–573 (2010). 771 25. Marlowe, F. W. Hunter-gatherers and human evolution. Evol. Anthropol. Issues, News, 772 Rev. 14, 54–67 (2005). 773 26. Burger, J. R. & Fristoe, T. S. Hunter-gatherer populations inform modern ecology. Proc. 774 Natl. Acad. Sci. 115, 1137–1139 (2018). 775 27. Hurtado, A. M. & Hill, K. R. Seasonality in a Foraging Society: Variation in Diet, Work 776 Effort, Fertility, and Sexual Division of Labor among the Hiwi of Venezuela. J. 777 Anthropol. Res. 46, 293–346 (1990). 778 28. Wilmsen, E. N. Studies in diet, nutrition, and fertility among a group of Kalahari 779 Bushmen in Botswana. Soc. Sci. Inf. 21, 95–125 (1982). 780 29. Lee, R. B. What Hunters Do for a Living, or, How to Make Out on Scarce Resources. in 781 Man the Hunter (eds. Lee, R. B. & DeVore, I.) 30–48 (Aldine de Gruyter, New York, 782 1968). 783 30. Hamilton, M. J., Milne, B. T., Walker, R. S. & Brown, J. H. Nonlinear scaling of space 784 use in human hunter-gatherers. Proc. Natl. Acad. Sci. 104, 4765–4769 (2007). 785 31. Messer, E. Anthropological Perspectives on Diet. Annu. Rev. Anthropol. 13, 205–249 786 (1984). 23 787 32. Testart, A. et al. The Significance of Food Storage Among Hunter-Gatherers: Residence 788 Patterns, Population Densities, and Social Inequalities [and Comments and Reply]. Curr. 789 Anthropol. 23, 523–537 (1982). 790 33. Winterhalder, B. Diet choice, risk, and food sharing in a stochastic environment. J. 791 Anthropol. Archaeol. 5, 369–392 (1986). 792 34. Kelly, R. L., Pelton, S. R. & Robinson, E. Studying sharing from the archaeological 793 record: problems and potential of scale. in Towards a Broader View of Hunter-Gatherer 794 Sharing (eds. Lavi, N. & Friesem, D. E.) (2019). 795 35. Joannes-Boyau, R. et al. Elemental signatures of Australopithecus africanus teeth reveal 796 seasonal dietary stress. Nature 572, 112–115 (2019). 797 36. Smits, S. A. et al. Seasonal cycling in the gut microbiome of the Hadza hunter-gatherers 798 of Tanzania. Science 357, 802–806 (2017). 799 37. Barnosky, A. D. Assessing the Causes of Late Pleistocene Extinctions on the Continents. 800 Science 306, 70–75 (2004). 801 38. Henrich, J. Demography and Cultural Evolution: How Adaptive Cultural Processes Can 802 Produce Maladaptive Losses—The Tasmanian Case. Am. Antiq. 69, 197–214 (2004). 803 39. Powell, A., Shennan, S. & Thomas, M. G. Late Pleistocene Demography and the 804 Appearance of Modern Human Behavior. Science 324, 1298–1301 (2009). 805 40. D’Alpoim Guedes, J. A., Crabtree, S. A., Bocinsky, R. K. & Kohler, T. A. Twenty-first 806 century approaches to ancient problems: Climate and society. Proc. Natl. Acad. Sci. U. S. 807 A. 113, 14483–14491 (2016). 808 41. Cegielski, W. H. & Rogers, J. D. Rethinking the role of Agent-Based Modeling in 809 archaeology. J. Anthropol. Archaeol. 41, 283–298 (2016). 810 42. Axtell, R. L. et al. Population growth and collapse in a multiagent model of the Kayenta 811 Anasazi in Long House Valley. Proc. Natl. Acad. Sci. 99, 7275–7279 (2002). 812 43. Hayden, B. Research and Development in the Stone Age: Technological Transitions 813 among Hunter-Gatherers. Curr. Anthropol. 22, 519–548 (1981). 814 44. Itkonen, T. I. Suomen lappalaiset vuoteen 1945. Ensimmäinen osa. (WSOY, Porvoo, 815 Finland, 1848). 816 45. Kirby, K. R. et al. D-PLACE: A Global Database of Cultural, Linguistic and 817 Environmental Diversity. PLoS One 11, e0158391 (2016). 818 46. NTSG (Numerical Terradynamic Simulation Group). MODIS NPP (MOD17A3). 819 Available at: http://files.ntsg.umt.edu/data/NTSG_Products/MOD17/. (Accessed: 12th 820 March 2015) 821 47. Defries, R.S., M.C. Hansen, F.G. Hall, G.J. Collatz, B.W. Meeson, S.O. Los, E.Brown De 822 Colstoun, and D.R. Landis. 2009. ISLSCP II Continuous Fields of Vegetation Cover, 823 1992-1993. ORNL DAAC, Oak Ridge, Tennessee, USA. 824 https://doi.org/10.3334/ORNLDAAC/931. 24 825 48. Bodesheim, P., Jung, M., Gans, F., Mahecha, M. D. & Reichstein, M. Upscaled diurnal 826 cycles of land–atmosphere fluxes: a new global half-hourly data product. Earth Syst. Sci. 827 Data 10, 1327–1365 (2018). 828 49. Šímová, I. & Storch, D. The enigma of terrestrial primary productivity: measurements, 829 models, scales and the diversity-productivity relationship. Ecography (Cop.). 40, 239–252 830 (2017). 831 50. Bontemps, S. et al. Consistent Global Land Cover Maps For Climate Modelling 832 Communities: Current Achievements Of The ESA’ Land Cover CCI. ESA Living Planet 833 Symp. (2013). 834 51. Bliege Bird, R. & Bird, D. W. Why Women Hunt - Risk and Contemporary Foraging in a 835 Western Desert Aboriginal Community. Curr. Anthropol. 49, 655–693 (2008). 836 52. Bliege Bird, R., Codding, B. F. & Bird, D. W. What Explains Differences in Men’s and 837 Women’s Production? Hum. Nat. 20, 105–129 (2009). 838 53. Reyes-García, V., Díaz-Reviriego, I., Duda, R., Fernández-Llamazares, Á. & Gallois, S. 839 “Hunting Otherwise”-Women’s Hunting in Two Contemporary Forager-Horticulturalist 840 Societies. Hum. Nat. 31, 203–221 (2020). 841 54. Krinner, G. et al. A dynamic global vegetation model for studies of the coupled 842 atmosphere-biosphere system. Global Biogeochem. Cycles 19, (2005). 843 55. Hamilton, M. J., Lobo, J., Rupley, E., Youn, H. & West, G. B. The ecological and 844 evolutionary energetics of hunter gatherer residential mobility. Evol. Anthropol. Issues, 845 News, Rev. 25, 124–132 (2016). 846 56. Abrams, P. A. & Ginzburg, L. R. The nature of predation: prey dependent, ratio dependent 847 or neither? Trends Ecol. Evol. 15, 337–341 (2000). 848 57. Winterhalder, B., Baillargeon, W., Cappelletto, F., Randolph Daniel, I. & Prescott, C. The 849 population ecology of hunter-gatherers and their prey. J. Anthropol. Archaeol. 7, 289–328 850 (1988). 851 58. Illius, A. W. & O’Connor, T. G. Resource heterogeneity and ungulate population 852 dynamics. Oikos 89, 283–294 (2000). 853 59. Golley, F. B. Energy Values of Ecological Materials. Ecology 42, 581–584 (1961). 854 60. Herbers, J. M. Time resources and laziness in animals. Oecologia 49, 252–262 (1981). 855 61. Raichlen, D. A. et al. Sitting, squatting, and the evolutionary biology of human inactivity. 856 Proc. Natl. Acad. Sci. 201911868 (2020). doi:10.1073/pnas.1911868117 857 62. Abrams, H., Jr. 1987 The Preference for Animal Protein and Fat: A Cross-Cultural 858 Survey. In Food and Evolution, edited by M. Harris and E. Ross, pp. 207–23. 859 Philadelphia: Temple University Press. 860 63. Hanya, G. & Aiba, S. Fruit fall in tropical and temperate forests: implications for 861 frugivore diversity. Ecol. Res. 25, 1081–1090 (2010). 862 64. Gherardi, L. A. & Sala, O. E. Global patterns and climatic controls of belowground net 863 carbon fixation. Proc. Natl. Acad. Sci. 202006715 (2020). doi:10.1073/pnas.2006715117 25 864 65. van Zonneveld, M. et al. Human diets drive range expansion of megafauna-dispersed fruit 865 species. Proc. Natl. Acad. Sci. 115, 3326–3331 (2018). 866 66. Max Roser, Hannah Ritchie and Esteban Ortiz-Ospina (2020) - ‘World Population 867 Growth’. Published online at OurWorldInData.org. Retrieved from: 868 “https://ourworldindata.org/world-population-growth” [Online Resource]. 869 67. Pontzer, H. et al. Metabolic acceleration and the evolution of human brain size and life 870 history. Nature 533, 390 (2016). 871 68. Viovy, N. CRUNCEP Version 7 - Atmospheric Forcing Data for the Community Land 872 Model. (2018). doi:10.5065/PZ8F-F017 873 69. Sobol′, I. . Global sensitivity indices for nonlinear mathematical models and their Monte 874 Carlo estimates. Math. Comput. Simul. 55, 271–280 (2001). 875 876 877 26 Population density (log10 ind. km-2) a NPP (gC m-2 yr-1) N=357, r=0.54 b log10(ΦNPP ) Figure 1 Energy intake FORGE (FORager dynamics by hunting in Global Ecosystem model) Time allocation Feedback Birth rate Birth rate Birth rate Mortality rate Mortality rate Mortality rate Browsers Grazers Energy expenditure Energy Excreta Energy intake (carbon) intake Energy (carbon) Excreta by gathering intake energy, water, Forests Fire Bare carbon cycling ground Grassland Multi-layer soil hydrology ORCHIDEE terrestrial and thermodynamics biosphere model Figure 2 a b density (log10 ind. km-2) Modelled population model obs Population density (ind./100km2) NPP (gC m-2 yr-1) c Foraging time d Gathering time e Hunting time hour day-1 ind.-1 Figure 3 a c Observation N=340, r=−0.60 ΦNPP Modelled meat fraction of the diet Meat fraction of the diet b d Model model obs ΦNPP Modelled log10(ΦNPP) Meat fraction of the diet Figure 4 Gathering-dominated regions 95.0% Gathering Plant into the soil pool Hunting Meat Spring 100% NPP (685 g C m-2 yr-1) 95.8% into the 0.76% 4.96% litter pool consumed by respired or excreted herbivores by herbivores Summer 4.2% 0.0022% consumed by consumed by hunter-gatherers Autumn herbivores 0.0079% Winter consumed by hunter-gatherers NPP (g C m-2 day-1) Food density (g C m-2) Hunting-dominated regions 93.0% into the soil pool 100% NPP (346 g C m-2 yr-1) Spring 94.0% 0.97% into the consumed by 6.97% litter pool herbivores respired or excreted by herbivores Summer 6.0% consumed by 0.0021% Autumn herbivores consumed by hunter-gatherers Winter 0.0017% consumed by hunter-gatherers NPP (g C m-2 day-1) Food density (g C m-2) Figure 5