Introduction

The Middle Pleistocene was a critical period for human evolution in Europe, marked by cyclic climate instability. Since 640 ka onwards shifts in the climate occurring every 100 ka caused a succession of glaciations and interglacials that are usually referred to the Marine Isotope Stages (MIS) and substages1. The four stages correlated with the 560 ka to 360 ka interval are MIS14 to MIS11. The MIS12 (Anglian or Elsterian glaciation) is considered a major turning point for the human occupation of Europe, marked by an increased abundance of archaeological data across Eurasia after this period2. Cultural complexity increased in Western Europe from MIS12 onwards, with the onset of Levallois, the spread of bifacial technology, and more complex resource management2,3. Moreover, the palaeontological evidence and many palaeogenetic analyses suggest that the earliest Neanderthal features appeared across Western Europe around 600–450 ka (roughly MIS15–MIS12), and the evolutionary process that gave rise to Neanderthals seems to have been complex4,5,6,7,8,9. Most models dealing with cultural and biological changes in Middle Pleistocene Europe include the effect of repeated shifts in the latitudinal distribution of humans, promoted by the cyclic climate oscillations characteristic of this period. It is generally assumed in these models that the distribution range of the European hominins expanded northwards during the interglacial stages and collapsed towards the Iberian, Apennine (or Italian), and Balkan Peninsulas during the glacial stages1,8. However, a precise delimitation of the range occupied by humans during each glacial and interglacial stage is currently impossible because of the relative scarcity of sites from this period. Arguably, the range expansions and contractions affected not only the spatial distribution of hominins but also the sizes of their census and effective populations. It is debatable whether increments in population size affected cultural changes in the Palaeolithic10,11,12, but oscillations in population size likely influenced biological evolution by modifying the frequency of some alleles and favouring the spread of certain genotypes by random genetic drift13.

Several different approaches have been developed to estimate the population density of Palaeolithic groups. These approaches include estimates derived from comparing the relative frequencies of different taxa between fossil assemblages and recent ecosystems14, extrapolation of ethnographic data15 combined with the observed frequency of archaeological sites16,17,18 or with the estimated size of the inhabitable area19, mathematical models simulating resource availability and competition in past environments20, and climate envelope models21. In addition to estimates of ancestral population size based on genetic data22, models producing continuous predictions of population density over a continent have proven a feasible approach for estimating past population sizes21. However, estimations of population size specific to the MIS14–MIS11 period in Europe are not available.

This study aimed to make two key contributions to the debate on the dynamics of cultural and biological evolution in Middle Pleistocene Europe. First, we present estimations of the potential human distribution in Western Europe during 11 cold and warm intervals of the MIS14–MIS11 period23 based on the lettering scheme for marine isotope substages24. Second, we provide human population size estimations that would have been sustained by Western European ecosystems during these 11 intervals. Our focus was on Western Europe because it has the most complete palaeontological and archaeological records for the Middle Pleistocene and, consequently, most models of cultural and biological evolution were developed specifically for this area2,8,15,25.

We applied an ecological niche modelling approach to estimate the potential distribution of humans for each interval based on a dataset of 68 archaeological assemblages from Western Europe, confidentially dated to this period (see “Methods” and Fig. 1). Palaeoclimate data for each interval of the MIS14–MIS11 period were obtained from the Oscillayers palaeoclimatic database26,27. Following other studies that modelled human distributions from archaeological data28,29, we assumed the Eltonian hypothesis30, which considers that the distribution of a species, hominins in our case, is determined by abiotic factors and minimally affected by biotic interactions. Thus, we selected the six most informative predictors from 19 bioclimatic variables and elevation (see “Methods”).

Figure 1
figure 1

Overview of the methodological approach for estimating total sustainable human populations in Western Europe from MIS 14 to MIS11. (a) Population density of recent hunter-gatherer groups and net primary productivity (NPP), both obtained from published sources (see “Methods”), were combined to obtain predictive equations of population density from NPP. 19 bioclimatic variables obtained from the Oscillayers dataset (see “Methods”) and elevation obtained from a digital elevation model were used as predictors in a species distribution model (SDM) to obtain a map of the suitable area for humans during each interval (see “Methods”). Presence data for the SDM were obtained from a dataset of archaeological sites dated to the period 560 ka to 360 ka, NPP maps for the interval were obtained from mean annual temperature (BIO1) and annual precipitation (BIO12) by applying the Miami model (see “Methods”). (b) The following procedure was repeated for each interval. Maximum and minimum population density maps were obtained by applying the predictive equations to the NPP map. Eventually, the total sustainable population in the interval was estimated by multiplying the average population density inside the predicted distribution range of humans by the size of the area (see “Methods”). The maps were created in QGIS 3.22 (https://www.qgis.org).

Net primary productivity (NPP) is a major determinant of hunter-gatherer population density, especially in the Holarctic31. Thus, we used a dataset of recent hunter-gatherer populations32 from Holarctic zones and NPP data31 to obtain a predictive equation for the maximum sustainable hunter-gatherer population density. It should be noted that, although our approach is a simple projection of the relationship observed between ethnographic hunter-gatherer population density and NPP, there are strong theoretical reasons to assume that there is a causal relationship between environmental productivity and hunter-gatherer population density31,33. The population density at which a species may occur in a particular environment is limited by the availability of trophic resources, which is directly dependent on NPP34. Since it is a well-known ecological rule that secondary consumers in natural communities have lower population densities than primary consumers35,36, a different relationship between population density and NPP in hunter-gatherer groups relying on different procurement strategies may be expected. Recent hunter-gatherer groups may be classified as ‘hunters’, ‘gatherers’ or ‘fishers’ according to their main procurement strategy32, although most of them use all three types of resources to some degree. Although the oldest evidence of fish consumption dates back to 1.95 Ma.37, foragers would have required more advanced technology than was available during the Lower Palaeolithic38 to rely on fishing as the main procurement strategy. Indeed, all archaeological assemblages dated to the Lower Palaeolithic that have faunal remains have provided overwhelming evidence for mammal consumption. In contrast, the evidence of aquatic resource consumption is very scarce for sites older than 160 ka39. Furthermore, well-preserved wooden hunting spears dated to around 300 ka were found at Schöningen40, but there is no direct or indirect evidence of nets or harpoons until the Upper Palaeolithic. Moreover, aquatic resources seem to have substituted for meat for contemporary hunter-gatherers who mainly rely on aquatic resources41. Thus, we excluded ‘fishers’ from our model to estimate population density in the Lower Palaeolithic and explored the relationships between NPP and density for ‘hunters’ and ‘gatherers’ separately (see “Methods”).

We estimated the sustainable population of Western Europe during each interval by multiplying the average hunter-gatherer sustainable population density in the area with suitable environmental conditions for humans by the size of the area (Fig. 1b).

Results

Relationship between NPP and population density

As already noted by other authors31, there is a significant relationship between NPP and population density for the entire sample of hunter-gatherer groups (Spearman’s rho = 0.68; p < 0.05, N = 215). Interestingly, the relationship is markedly different based on the main subsistence strategy (Table 1, Fig. 2). Excluding horseback hunters from the dataset did not significantly change the slope of the relationship (F = 1.099, p = 0.7291; t = − 1.061 p = 0.2911, N = 101, see Supplementary Figure S1).

Table 1 Simple linear regression of log-transformed hunter-gatherer population density (D) on log-transformed net primary productivity (LOGNPP). Regression equations were fitted to the whole sample (average) and the upper and lower limit of the distribution for each procurement strategy.
Figure 2
figure 2

Relationship between net primary production (NPP) and hunter-gatherer population density (D). Green squares represent groups classified as ‘gatherers’ and red circles correspond to ‘hunters’. The solid lines are the regression equations fitted to the total sample of ‘gatherers’ (green) and the total sample of ‘hunters’ (red). Dashed lines are the regression equations to the upper and lower limit of the distribution of ‘gatherers’ (green) and ‘hunters’ (red), respectively.

The slopes of the lines show that ‘gatherers’ were able to sustain higher population densities than ‘hunters’ at a given amount of NPP. This aligns with the observation that ‘hunters’ require larger home ranges than “gatherers” 33.

Ecological Niche model

Comparing different model candidates highlighted the necessity of calibrating the species distribution model (SDM) for individual species42,43,44. Our results (Fig. 3) showed that more complex feature classes led to considerably better performance in validation metrics, such as the increment in the Akaike information criterion for small sample sizes (∆AICc) and the area under the receiver-operator curve (AUC). However, we interpreted this as an effect of overfitting, because the complex models showed much higher omission rates in the testing subsample. Therefore, we prefer a compromise model that delivers reasonable values for all indicators considered. Furthermore, we could not find ‘peak performance’ at the default regression multiplier (RM) of 1.0, but rather observed a trend indicating that lower RMs produced better results than higher RMs. This observation contradicts other studies43,44,45 that recommend higher RMs, which impose higher penalties to include covariates, and thus, simplify the model. However, such studies used larger and unreduced covariate sets that did not undergo a preselection process (see “Methods”), so collinear and correlated covariates were likely penalised through the higher RM. We interpreted this as an advantage of our statistical preprocessing, which retained only covariates with high individual information content, and thus, were able to inform more complex models with lower penalties.

Figure 3
figure 3

The human ecological niche throughout MIS 11–14 according to the selected SDM. (a) Comparison of densities at confirmed presence points (n = 644; red) with a random sample (n = 10,000; grey) of the overall environment. (b) Depictions of the human niche (red) within the environment (grey) in a covariate space with the 3 most important covariates. (c) Ranking of covariate importance based on the metrics permutation importance (blue) and model contribution (yellow).

The most informative variables for niche estimation were annual precipitation (BIO12), precipitation seasonality (BIO15), and mean annual temperature (BIO1), while mean diurnal temperature range (BIO2), elevation, and temperature annual range (BIO7) played minor roles (Fig. 3c). Thereafter, the major limiting factors for the human niche were extremely dry (BIO12 < 300 mm) and wet (BIO12 > 1200 mm) conditions, strong precipitation seasonality (BIO15 > 65), and mean annual temperatures below 0 °C (Fig. 3a and b).

The niche model showed great stability in the extent and distribution of the permanent core areas suitable for hominins from MIS14 to MIS11 and a temporally varying periphery (Supplementary Figure S2). The most variable areas corresponded to the North Sea, the Jutland Peninsula, and adjacent territories. Moreover, the region east of the Oder River was only suitable for humans during MIS11c (Fig. 4). This variability was determined by the unsuitability of the northern region during the glacial stages, especially during MIS12. The niche model predicts a continuous suitable area connecting the southern peninsulas to Britain through France and West Germany for most of the time. Britain’s connection with the continent was only interrupted during the MIS11c and MIS13a. The areas predicted to be non-suitable throughout the entire MIS14 to MIS11 period were the high-altitude regions of the Pyrenees and the Alps; Brittany; two small spots in the Ardennes and central France; and a large part of the Iberian Peninsula, including the southern half and the northwest corner. The exclusion of a large area of the Iberian Peninsula might be explained by the unique climatic conditions of this region and the absence of sites confidentially dated to this period.

Figure 4
figure 4

Variation in the area suitable for humans during the MIS14 to MIS11 period. The maps were created in QGIS 3.22 (https://www.qgis.org).

Total population from MIS14 to MIS11 in Western Europe

The estimated sustainable population density in Europe is distributed following a consistent pattern throughout the MIS14 to MIS11 periods. A marked West to East gradient is evident, with high values along the Atlantic coast that diminish towards the interior of the continent. This longitudinal gradient is especially relevant in the Iberian Peninsula, where the sustainable population density also shows high variation over time along the Atlantic coast compared with the interior and eastern areas (Fig. 5).

Figure 5
figure 5

Population density in Western Europe from MIS14 to MIS11. The maps were created in QGIS 3.22 (https://www.qgis.org). in the areas considered suitable for humans based on the SDM. The maps illustrate the lower limit of the estimated interval for population density.

The estimated maximum total population along Western Europe shows little variation throughout the MIS14 to MIS11 interval. The maximum value obtained for the MIS11c interval was between 25,000–170,000 individuals while the minimum value for MIS12b was between 13,000–85,000 individuals. This estimate represents a 50% decrease in the potential population at the hardest glacial interval in comparison with the most favourable interglacial interval. Interestingly, the variations in the potential total population of Western Europe were driven by changes in the extent of the area suitable for humans, not by variations in the maximum sustainable density at the local scale (Fig. 6). Indeed, the average population density was constant throughout this period, while the periphery of the suitable area showed marked variation.

Figure 6
figure 6

Changes in average population density (D), total population and suitable area in Western Europe during the MIS14-MIS11 period. The LR04 stack of marine benthic foraminiferal δ18O46 is shown as a proxy for climate variation, indicating the substages defined for the MIS14–MIS11 period 24.

Discussion

It is generally assumed that climate cyclicity severely impacted human distribution in Western Europe during the Middle Pleistocene by restricting the suitable area for humans to the southernmost region of the continent1,3,8,25,47. However, the results of the ecological niche model presented here suggest a different scenario. Although most of Britain and large areas north of parallel 50° N were unsuitable for humans during the coldest periods of the MIS14 to MIS11 interval, a core area with favourable conditions for human occupation spread continuously from northern Iberia, Italy, and most of France. Thus, the area suitable for humans during cold periods would have been much larger than normally assumed and would not have been restricted to the Mediterranean peninsulas. Consequently, the populations inhabiting the northwestern area were not entirely disconnected from southern populations as a result of climate shifts, as proposed by some models of human evolution in Europe1,8,25. The results of the SDM suggest the existence of a continuous environmental corridor to the east of the Pyrenees during both the cold and warm phases. This corridor would have connected the Iberian and North Western human populations throughout the entire MIS14 to MIS 11 period. In contrast, the Apennine peninsula was disconnected from the rest of the core area during the coldest stages (the entire MIS12 and MIS14 glacial stages). Consequently, human populations in the Apennine peninsula would have been isolated from other populations in Europe during these cold stages, which might have favoured their genetic and cultural divergence from other European populations.

Most scholars agree that Britain was not populated by humans during MIS 12 and MIS 148,15,48. However, the SDM predicts the presence of suitable areas for humans in southern Britain during cold phases. Lithic artefacts have been retrieved from sediments dated to the MIS12 and MIS14 glacial stages in the Bytham River area, but they are usually regarded as reworked items dating back to the former interglacial stage49. If this interpretation is correct, and Britain was depopulated during the glacial stages, the SDM might have overestimated the area suitable for humans during those periods. It must be noted, however, that an SDM produces a map of the areas potentially suitable for the target species based on the factors included in the model. In our case, these factors were climate and elevation. Therefore, that area represents an approximation of the fundamental niche50 of the Middle Pleistocene European hominins. Although the fundamental niche is generally restricted to biological adaptations, we include here also the potential cultural adaptations of the Middle Pleistocene humans in our use of the term. The realised niche is determined by the overlap of the fundamental niche and the region with an appropriate combination of biotic factors (food resources, competitors, predators, etc.) for the target species50. Finally, the region in which a species is actually present is the geographic extent of the realized niche accessible for the species50. Thus, even if humans were able to cope with the climatic conditions of Britain during MIS14 and MIS1223, other factors, such as the lack of suitable food resources, may have excluded them from that area.

Another unexpected result of the SDM was the exclusion of the Atlantic Coast and the southern half of Iberia from the geographic projection of the fundamental niche of the Middle Pleistocene European hominins for most of the MIS14 to MIS11 period. The Atlantic Coast is significantly more humid, and the southern half is considerably drier and warmer than the rest of the peninsula. There are no other similar climate areas in Western Europe, and there are no sites in that area with a chronology precise enough to be included in the model. Thus, this result should be considered cautiously, as it could change with the discovery of a single site from southern or western Iberia dated to this period.

According to our results, the changes in the total human population that could be sustainable by the environment in Western Europe were largely determined by the variation in the size of the suitable area during the MIS14 to MIS11 period (Fig. 6). The sustainable population density was constant over time and exhibited little spatial variation. Consequently, the effect of changes in population density on the size of the total sustainable population was negligible. The small temporal and spatial variation in population density was a direct consequence of the slight variation in NPP, in time and space, within the area having suitable conditions for humans. It may be argued that our methodology overestimated the amount of resources available to humans, especially in harsh environments, because the Miami model estimates total net primary productivity and only a small amount of that productivity may be directly available to humans. However, our method does not make assumptions about the proportion of NPP actually used by humans. The estimated values of population density are based on the empirical relationship between the population density of recent hunter-gatherers in a given area and the NPP estimated for that area by the Miami model. The proportion of total NPP that can be used by humans varies across different biomes32. Thus, it may be expected that the proportion of NPP available to humans in Western Europe was different during glacial and interglacial periods, but this fact does not affect the estimation of population density. The database used to obtain the predictive equations included recent hunter-gatherer populations from a wide sample of biomes, ranging from desert and prairie to taiga and tundra. Thus, it includes the main biomes occurring in Western Europe during the cold and warm phases of the MIS14 to MIS11 periods51,52.

Our approach assumes that the population density of Middle Pleistocene hominins may be inferred from ethnographic hunter-gatherer population densities. This assumption may be negatively affected by the biological, cultural and technological differences that exist between the recent hunter-gatherer populations used to derive the empirical relationship between NPP and population density, and the Lower Palaeolithic hunter-gatherers. It has been proposed that Neanderthals lived in smaller groups than recent hunter-gatherers based on genetic data22. Differences in social organization, derived from differences in their cognitive capacities, might force Middle Pleistocene humans to live also in small groups and small social networks. This would have an impact on their population density. It has been proposed that group size is related to brain size in primates53, but see also54, and the brain size of the Middle Pleistocene European hominins was similar to, but smaller than, the brain size of Homo sapiens55. Moreover, recent hunter-gatherers use sophisticated hunting weaponry and a set of food conservation techniques not available to Lower Palaeolithic hunter-gatherers. These differences likely make recent hunter-gatherers more efficient than Lower Palaeolithic hunter-gatherers in extracting and processing resources, allowing them to co-opt a larger proportion of NPP and achieve higher population densities than Lower Palaeolithic hunter-gatherers in an equivalent environment. This is a major factor that may have inflated the estimated population density values for the Middle Pleistocene populations. However, the large similarities in body size, diet, and physiology between Middle Pleistocene hominins and H. sapiens support the validity of the aforementioned assumption. Therefore, taking the estimated values of population density with caution, only the lower limit of the estimated interval for population density will be considered in the following discussion.

Following this conservative approach, the total sustainable population in Western Europe during the MIS14–MIS11 interval would have oscillated between 13,000 and 25,000 individuals. Models of human evolution in Europe during the Middle Pleistocene often estimate a census population five to ten times lower8, giving an effective breeding population as low as 600–1,000 individuals in full glacial times and 1,200–2,000 individuals during warm periods8. There is strong theoretical and empirical evidence suggesting that the optimal size for a hunter-gatherer group or band is 25 individuals41, although the validity of the multi-scalar organization pattern commonly assumed for hunter-gatherer societies has recently been questioned54. Assuming an average size of 25 individuals for a Lower Palaeolithic hunter-gatherer band, it has been estimated that there may have been 50–100 such bands in Europe during full glacial periods and 120–200 bands during interstadials8. According to the widespread assumption that the European population was confined to refugia in the Iberian, Apennine, and Balkan Peninsula during the glacial stages1,3,8,25, only 20–30 such groups would have been isolated within each peninsula. In this scenario, the probability of interbreeding between groups of the same peninsula would be very low, and the genetic flow between populations in different peninsulas would have been impossible for several millennia. Moreover, assuming that the total census population was 1,500–2,000 individuals, with only 40% being of reproductive age8; that interbreeding between populations from different peninsulas was not possible; and, for the sake of simplicity, that the total population was equally distributed among the three areas; the effective breeding size of each isolated population would have been only 200–333 individuals. Under these assumptions, the census population would be 500–833 for each peninsula, which is much lower than the minimum viable population size for any organism, which exceeds a few thousand individuals56. In contrast, under the new scenario described here, the number of bands in Western Europe could have been as many as 520 in full glacial periods and up to 1,000 during warm periods. Moreover, the Iberian bands were connected to the north of the Pyrenees during MIS12 and MIS14, allowing genetic flow between them and the bands in the north. Only the population from the Apennine peninsula would have been isolated from the others; however, according to the minimum sustainable population density estimated here for that area during MIS14 and MIS12, this population could have exceeded 2,000 individuals. A population of that size isolated in the Apennine peninsula over several millennia could be viable, although it would also have a high risk of extinction. In contrast, the main metapopulation inhabiting France and Iberia could attain a size far above the threshold for a viable population during the entire MIS14 to MIS11 period.

The estimated values of the total population reported here represent the maximum number of individuals that the Western European ecosystems could sustain during the MIS14 to MIS11 periods. These estimations represent an ecological baseline from which to think about Middle Palaeolithic population history in Western Europe. Other factors not included in our approach, like the cultural differences among groups, biological singularities of the Middle Pleistocene hominins, diseases or stochastic events, likely also played a role in determining the actual size of the Middle Pleistocene populations. Therefore, our results should be considered the maximum thresholds, and it must be acknowledged that the actual population size may have been lower. However, these numbers, combined with the geographic projection of the fundamental niche of Middle Pleistocene humans produced by the SDM, set up a new theoretical scenario for models and hypotheses aimed at explaining cultural and biological changes in Western Europe during the critical period of human evolution known as the Middle Pleistocene.

Methods

Occurrence data

This study used a published dataset of archaeological assemblages correlated with the interval from MIS14 to MIS1123 compiled from published literature. The dataset included archaeological assemblages from the European continent and the British Islands west of the meridian 14oE, along with the entire Apennine peninsula. No latitudinal limits were established within the European continent for sites to be included in the dataset. Each archaeological assemblage was correlated to the marine isotope curve24 according to information provided by the original sources on the biostratigraphy, chronostratigraphy and radiometric dates for the site23. The dataset includes 68 archaeological assemblages from 46 sites, but because of relatively imprecise dating methods available for this period, only 33 assemblages correlated to an interval of the marine isotope curve (Supplementary Table S1). The other 35 assemblages could not be correlated with certainty to an interval within a stage23.

Palaeoenvironmental data

Palaeoclimate data for the MIS14–MIS11 period were obtained from the Oscillayers palaeoclimatic database26,27. Oscillayers includes interpolated estimates for 19 bioclimatic variables (BIO1 to BIO19) at high spatial (2.5 arc-min) and temporal (10 ka time periods) resolution and covers the period from 5.4 million years ago to the present26. The set of 19 bioclimatic variables provided in Oscillayers is the same as that defined in the BIOCLIM model57: these variables are commonly used in ecological niche modelling analyses58. Because the duration of intervals varied from 10 to 30 ka, it was necessary in some cases to aggregate several Oscillayers time slices (Supplementary Table S2). When an interval overlapped the time period represented by two or more Oscillayer time slices, the values of the bioclimatic variables were obtained as the average 23. A map of NPP in grams of dry matter/m2·yr was estimated for each interval using the Miami model59 with the mean annual temperature (BIO1) and total annual precipitation (BIO12) obtained from the Oscillayers dataset.

All geodatasets in this study were transformed to the Lambert azimuth equal-area projection based on the ETRS89 ellipsoid. With a projection centre of Lon =  + 52° and Lat =  + 10 o, as defined in EPSG:3035, this allows high accuracy area measurements.

Ecological niche modelling

We chose to describe and infer the human niche across diverse climate regions during multiple warm and cold phases. For this, we selected statistically sound covariates to describe the environment and generated a sample accounting for the spatial surroundings of the sites. We performed computational tests on species SDMs using these data to assess the range of hominin habitable zones.

Regarding environmental factors, we started with 19 bioclimatic variables from the Oscillayer dataset and a digital elevation model. Elevation is usually deprecated as a proxy describing the landscape60. However, it was included here to evaluate if adaptation to height might play a role in the distribution of the target hominin populations. More sophisticated topographic variables derived from high-resolution DEMs are often involved in studies targeted at site prediction, where they exert strong predictive power. The distribution of early archaeological sites is often related to certain types of geology and topography (e.g. limestone caves, gorges, undercuts) due to site preservation factors. However, the aim here is not the prediction of sites but producing an approximation to the human niche. Therefore, we excluded the complex and biased factor ‘landscape’ from our study. Because the Maximum Entropy (MaxEnt) algorithm provides functions for expanding the covariates (i.e. the input environmental variables) into complex feature classes, it is recommended to reduce the number of variables in advance to prevent undesirable effects caused by correlation and multicollinearity61,62. A common approach to reduce dimensionality is to project the variables onto principal components63,64. However, this comes at the price of a reduced SDM interpretability or highly noncausal relationships65. Instead, we selected meaningful covariates through thresholds. Therefore, we calculated a pairwise correlation matrix of all variables with the R package “Hmisc” v.4.566 and systematically excluded those that were correlated with others using a Pearson or Spearman correlation coefficient of |r|> 0.9. Where multiple correlated variables existed, we retained the more general variable (e.g. “Annual Mean Temperature” over “Mean Temperature of Coldest Quarter”). We further assessed the issue of multicollinearity between variables by calculating R2, Tolerance, and Variance Inflation Factor (VIF) with the package “fuzzySim” v.3.067. These indicators describe the linear dependence of one variable on multiple other variables (Supplementary Table S3 and Supplementary Figure S3). Since we found concerning amounts of collinearity, we iteratively reduced the variables with the highest VIF until only those variables with a threshold VIF < 5 remained. Recommended thresholds lie at VIF < 10, VIF < 5, or VIF < 368,69. Through this selection process, the candidate variables elevation, BIO1, BIO2, BIO7, BIO12 and BIO15 were chosen as suitable covariates for further modelling.

The locations of archaeological finds are often related to specific geologic and taphonomic conditions and do not necessarily reflect the environmental variability of the landscape that allowed early humans to survive. However, these surroundings are elementary for hunting and gathering lifestyles and need to be addressed in an SDM. Therefore, we expanded the occurrence site locations by including their geographical surroundings within a buffer radius of 10 km from the sample. This radius corresponds to maximal one-way distances of daily foraging trips of modern hunter-gatherers41 and resulted in a total sample size of 664 points around 30 sites. Where sites were occupied multiple times and assemblages lay in different substages, these were treated as individual sites with their respective environmental conditions. We accounted for spatial autocorrelation by creating a grid that allowed only one sample per raster cell. Furthermore, only onshore cells were considered, based on the palaeo sea-level estimates within the Oscillayer dataset, which rely on the sea-level reconstructions and the ETOPO1 elevation model27,70.

Because our archaeological dataset is based on occurrences and lacks confirmed absence data, we chose the MaxEnt algorithm, which is considered a powerful tool for SDM with presence and background data71,72,73. Our analyses were conducted with the software ENMeval v.2.0.142, which relies on MaxEnt v.3.4.474 and is available in the statistical software R v.4.0.575. The SDM was calibrated by performing tests with the following alternate configurations:

  • The feature space was enhanced by applying six functions (L, LQ, LQH, LQP, LQHP, and LQHPT with L = linear, Q = quadratic, H = hinge, P = product, T = threshold) to the initial covariates.

  • We tested eight regularisation multipliers (0.2, 0.4, 0.6, 0.8, 1.5, 2, 3, and 4) in addition to the default value of 1.0.

  • Because MaxEnt is capable of penalizing covariates and because of our rigid preselection procedure, we did not perform additional tests with varying covariate sets.

Each of these experiments was tested with the same 10,000 background points, initial covariates and sample data. A sensitivity analysis showed that the distribution of the sampled values can be considered stable when the sample exceeds the size of n > 1,000 (Figure S4a). The models were run with 500 iterations and a k-fold subsampling strategy to ensure that sites and their surroundings were treated as a whole. In earlier attempts, we noticed that completely random subsets resulted in heavily over-fitted models; this was attributed to spatial autocorrelation from our sampling design that characterized single sites by multiple points.

We evaluated these models based on three indicators (Fig. 7), each computed with ENMeval: the Akaike information criterion for small sample sizes relative to the “best” performing model (∆AICc), the average Omission Rate of the 10%-percentile of presence points (OR.10), and the area under the receiver-operator-curve (AUC) of the test subsample76. The model performance results showed, that the feature set choice (indicated by colour) majorly affected the model quality. Particularly, those models that outperformed in one indicator performed poorly in others. For example, the feature set LQHPT performed best in AICc but worst in OR.10, whereas feature set L had the best performance in OR.10 but poor values regarding AICc and AUC. Therefore, we chose a compromise model, with a feature set LQP and a regularisation multiplier of 0.4. This model combined an ∆AICc = 624 (AICc = 10,504) lying in the midfield, a good OR.10 = 0.14, and a high AUC = 0.86.

Figure 7
figure 7

Performance of model candidates evaluated with the indicators Akaike Information Criterion for small samples (AICc), Omission Rate of the testing points at 10% training threshold and Area Under ROC Curve (AUC). The models varied in using feature classes (colour) and the regularisation multiplier (x-axis). The dashed vertical line indicates the default model with a default regularisation of 1.0.

We projected this model onto the environmental data of the 11 time periods and thus produced continuous surfaces with log transformations (1 = most suitable, 0 = least suitable), indicating suitability for human occurrences. These were further binarised to create a mask of the human habitat, based on a threshold of suitability = 0.118, corresponding to the 5%-percentile of predicted values at observed sites (see also28).

The 33 archaeological assemblages not correlated with certainty to an interval within a stage were used for further validation of the model (Supplementary Table S4).

Estimation of population density

A dataset of recent hunter-gatherer populations compiled by L.R. Binford32 and data on NPP31 were used to obtain a predictive equation of the maximum sustainable hunter-gatherer population density in individuals/100 km2 from NPP in grams of dry matter /m2·yr. The data on recent NPP are based on estimations from the Miami model21,59. This source was preferred to satellite-based values for two reasons. First, climate-based NPP estimates avoid the potential problems caused by recent human appropriations of NPP21, and second, we used the Miami model to estimate NPP in the past, making our estimates of NPP comparable through time. Latitude, longitude, population density, and percentage of fishing, hunting, and gathering in the diet of 215 Holarctic hunter-gatherer groups living above the parallel 30° N were obtained from Binford’s dataset32 (Supplementary Table S5). The population density was obtained by dividing the total population of the ethnic group by the total area occupied32. For the hunter-gatherer groups with additional population density estimates available41, we used the average of the two estimates. As highlighted by other authors21,32, the sample of recent hunter-gatherers in Binford’s dataset is geographically biased but is not biased concerning niche space or habitat. Hunter-gatherer societies considered dependent on trade with non-forager societies for their subsistence at the time of documentation32 were excluded. K-means clustering was used to classify the hunter-gatherer groups according to their procurement strategy using the percentage of hunting, fishing and gathering in their diet as input data (Supplementary Table S6). This analysis confirmed the validity of the classification of hunter-gatherer groups as ‘hunters’, ‘fishers’ and ‘gatherers’ provided by the original source32 based on their most frequent activity. However, because our results and the original source differed in the assignment of some hunter-gatherer groups to a procurement strategy group32, we retained the classification produced by k-means clustering in the subsequent analyses.

Population density and NPP were transformed to the logarithm of base 10 (LOGD and LOGNPP, respectively), and ordinary least square regression was used to fit a function to the relationship between LOGD and LOGNPP for the entire sample as well as for each separate procurement strategy group. One case was excluded from the analyses (Monachi-southern) after being identified as an outlier (out of the 95% prediction interval for the cases). Horseback hunting has been suggested to allow some hunter-gatherer groups to occur at higher densities than on-foot hunters at a given value of NPP21. Thus, we tested the effect of removing mounted hunters from the sample used to build the regression model. Regression equations were fitted to the upper and lower limits of the distributions of ‘Hunters’ and ‘Gatherers’33,77. The sample was divided into bins of 0.1 LOGNPP and the regression equations were fitted to the maximum and minimum population density for each bin. The lower and upper limit equations for ‘Hunters’ were applied to the NPP map for each time interval to obtain maps that estimated the maximum and minimum population densities at each cell. The maximum and minimum total populations in Western Europe for each time interval were estimated as:

$$P=\overline{D }\times A$$

where P is the total population, \(\overline{D }\) is the average population density (individuals/100 km2) within the range predicted for the substage and A is the size of that range in km2. Two P estimates were obtained for each time interval, the first from the maximum population density estimate and the second from the minimum population density estimate. K-means clustering and regression analyses were carried out using the StatSoft Statistica 13.