Timeframe of interest

The Last Glacial period (encompassing Marine isotope stages 4–2) is characterised by a global trend towards cooler and drier conditions, culminating with the Last Glacial Maximum (LGM), which correlates with Greenland Stadial 2.1 and marks the maximum volume of global ice and consequently, the lowest sea-stand [5–7]. Proxy-based climate data indicate that relatively cold and dry conditions prevailed in Europe at this time, accompanied by a decline in arboreal pollen and the expansion of the Eurasian steppe [8,9]. Various paleoclimate models produce comparable results [10–12]. The downward trend in temperature and precipitation rates that characterises the Last Glacial was punctuated by a series of abrupt, millenial-scale climate shifts, well-documented in the Greenland ice-core records [8,9]. These climate shifts may have played a role in the timing of modern human dispersals into Europe [13]. The onset of extreme conditions during the LGM is thought to have patterned the subsequent distribution of human populations within Europe [14–16].

The archaeological record indicates that human populations withdrew from Northern Europe and other previously inhabited regions, such as the centre of the Iberian Peninsula, during the LGM (Fig 1). Whether the populations living in these regions crashed or whether they migrated to glacial refugia is not currently known, but it seems likely that overall population numbers fell during this time interval [17]. Paradoxically, the record also shows that the LGM was a period of cultural innovation, marked by the regional fragmentation of material culture [18–20]. The advance of the ice-sheets and the retreat of human populations to southern, glacial refugia during the LGM hypothetically created a situation where smaller, regional populations became spatially isolated, producing new and divergent cultural patterns. This research allows us to address the question of the regionalisation of material culture in Western Europe during the LGM by modelling human population dynamics, including individual mobility, inter-regional connectivity and regional survival rates as a function of the structure of the environment, represented by a model of habitat suitability which acts as a proxy.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Spatial distribution of habitat suitability values derived in Burke et al. (2017).

Higher values of habitat suitability indicate cells that are more favourable for human habitation. Points represent known LGM period archaeological sites. Note that known glacial extents are masked. European border polygons republished from GADM.org under a CC BY license, with permission from GADM, original copyright 2018.

https://doi.org/10.1371/journal.pone.0217996.g001

Habitat suitability models

Underpinning this research is a model of habitat suitability for the LGM derived from the application of a machine-learning algorithm (Random Forest) to the spatial distribution of prehistoric sites and a suite of geographical and paleo-environmental variables [1]. Archaeologists use models such as this both predictively, e.g. to design archaeological surveys, and interpretively, as a means of identifying key predictors of human spatial behaviour. Climate variables such as precipitation rates and temperature are usually included in the modelling process because they regulate vegetation cover, which in turn determines net primary productivity (NPP), or available biomass. Geographical variables such as latitude, continentality, and elevation are included because they mitigate climate conditions and produce spatial and temporal variation in the distribution of resources key to human survival.

The advent of coupled global climate models combining atmospheric and oceanic circulation (GCMs) greatly enhanced our ability to produce accurate climate simulations of past and present conditions. The Paleoclimate Model Intercomparison Project (PMIP) was designed to test the ability of different GCMs to simulate climate under conditions of rapid or extreme climate change, focussing initially on the Last Glacial Maximum (centered at 21ka) and the mid-Holocene [21]. These paleoenvironmental simulations allowed scientists to compare and evaluate the relative performance of different model architectures, independently validating results using climate proxies derived from terrestrial and marine sediment records [16,21]. The IPSL-CM5 model used in this research compares very favourably with other GCMs participating in PMIP [12].

GCMs produce outputs on fine temporal but coarse spatial scales. Two methods exist for downscaling GCM outputs to a more appropriate scale for archaeological research: the use of nested regional dynamic models or statistical downscaling methods. The Stage 3 project, for example, which marks a watershed in the development of multidisciplinary research into human-environment interactions, used a conventional GCM coupled with a regional circulation model for Europe (REGCM2) to produce a meso-scale analysis (~60 km x 60 km) exploring the impact of climate change on Neanderthal and human populations during the last Glacial period [22,23]. The second approach, which is the one adopted by the Hominid Dispersal Research Group (HDRG), uses a statistical method to downscale the GCM outputs as a function of topography, continentality and temperature [24–26] producing simulated climate conditions at a fine spatial scale (15x15 km²). This approximates the scale of the “catchment” areas around archaeological sites, i.e., the area within which human activities took place on a daily basis, as suggested by Vita-Finzi and Higgs [27].

GCMs have been used elsewhere to produce spatial models of human spatial distribution during the LGM [28–33]. The various archaeological models that have been produced to date are broadly similar, which is understandable since they are based on essentially the same archaeological dataset (or subsets thereof) and comparable suites of candidate predictors. Predictors typically include elevation, slope, maximum and minimum temperatures and maximum and minimum precipitation rates. The devil is in the detail, however, and the habitat suitability model used in this research differs fundamentally from other models in its explicit use of climate variability (defined as deviation over a defined time-frame compared with long-term statistics) as a potential predictor of habitat suitability. Thus defined, climate variability is a measure of how predictable conditions are from year to year (the timeframe defined in our original study) with repercussions for the distribution of plant and animal resources. Numerical climate models are generally considered to underestimate interannual variability [34], at least under modern conditions, but still provide plausible physical scenarios with which to test the impact of variability on habitat suitability.

Hunter-gatherers use mobility as a key strategy for dealing with changes in the distribution of resources. Ethnographically known hunting and gathering groups maintain a fluid social structure and extended social networks that promote the gathering of information about resource availability and facilitate residential moves in the event of resource failure. Mobility has inherent costs that need to be balanced out against expected returns, however, and the structure of the environment, i.e., the pattern of distribution of both resources and resource risk, affects the degree of mobility in human systems [35,36]. Consequently, the HDRG team included several indices of climate variability in the model-building process in addition to a suite of conventional climate and topographic variables. The model-building methodology is described in full elsewhere [1]. In initial research at the scale of the Iberian Peninsula [16] and subsequent research on a larger European domain [1] climate variability consistently emerged as one of the predictors of archaeological site location. The habitat suitability model with the best performance (run 26) resulting from the 2017 publication is the one that is used in this research. According to this model, habitat suitability (HS) is determined by seven predictors: elevation, slope, variability in Spring precipitation rates, maximum monthly average temperature in Autumn, variability in monthly temperature calculated using the Standard Temperature Index, standard deviation in average monthly temperature in Winter, and the standard deviation in average monthly temperature in Autumn. Our working hypothesis is that Palaeolithic populations will have attempted to minimise ecological risk by avoiding regions where HS is relatively low.

Regional designations

Rather than model a homogenous landscape, we import a heterogeneous landscape based on the predictive model of habitat suitability (HS) (see above). This suitability landscape forms the landscape of the ABM. Instead of imposing regional population designations from previously published literature (such as Solutrean or Epigravettian for example), we use a spatial analysis procedure of the habitat suitability landscape itself to derive “hot-spots” of spatially clustered high HS-value cells. Specifically, we used Hot Spot Analysis Getis-Ord Gi* tool in ArcGIS 10.5 [39] with the distance threshold computed by the tool and then used all cells binned as +1 or above. We grouped the resulting clusters based on contiguous cells (Moore neighbourhood) using ArcGIS’s Region Group tool, and further grouped small isolated clusters of cells into larger groupings by hand digitizing large polygons with the borders running along a midline between each pair of clusters to produce six distinct clusters which we numbered sequentially (Fig 2A, see S1 Appendix for additional details). These clusters form the core areas of our “regions” and are the basis for attributing a genetic fingerprint to the agents (see below). The six “regions” (Fig 2B) are used as a heuristic device to track agent movements during the experiments (see below). Note that the model does not use the core area or region values to influence agent movement but only to track the movement and reproduction patterns of the metapopulation. Each 5 by 5 km cell in the ABM, therefore, has three attributes: habitat suitability (0–1), core area (1 to 6), and region (1 to 6).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

Example initial metapopulation of 7000 randomly distributed agents (coloured dots) overlaid on a) lighter coloured “core area” cells determined by the hotspot analysis, b) larger regions used for counting agent regional population sizes, and c) initial distribution of 7000 agents across the regions.

https://doi.org/10.1371/journal.pone.0217996.g002

The six core areas have similar HS values but are unequal in size, with the core area centered in France being far larger than the others (Table 1). The six regions have quite different cell frequencies across the range of habitat suitability values.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary data on the habitat suitability landscape per core area and region.

https://doi.org/10.1371/journal.pone.0217996.t001

Initial ABM model parameters

Population dynamics

We model the probability of agents disbanding or persisting, and the probability of creating new agents using life tables drawn from Pennington [43] and population growth rates derived from the ethnographic and archaeological literature. Pennington’s life tables contain the basic parameters that govern population growth rates. We use the birthrate calculated from these tables as a proxy for the probability that an agent may generate a new household.

To calculate the birthrate, we use Pennington’s [43] summary of the total number of births over a woman’s lifespan, i.e. the total fertility rate (TFR), for a variety of hunter-gatherer groups. Among the hunter-gatherer populations sampled, the TFR ranges from 2.6 to 8.0 births with a median around 4.3 (Table 7.4 in [43]). The hunter-gatherer groups in this sample live in a variety of environmental settings but the early Kutchin are a good proxy for LGM populations in Western Europe since they lived a traditional hunting and gathering lifestyle in subarctic conditions in the Canadian Yukon; their TFR is 4.4 (Ross 1981 cited in [43]). Although estimates of the age of first birth and last birth vary widely, the reproductive span of hunter-gatherer women extends from adolescence to about 40 years of age [43] (p. 175), therefore we assume an average reproductive span of 25 years. We took the Kutchin TFR of 4.4 and divided it by 25 to derive an annual birthrate of 17.6%.

The model runs with a time step representing one month. Over the course of each run, therefore, agents generate new households at a specified probability equal to birthrate divided by 12. This new birthrate parameter is the probability that agents will found a new household with the nearest other agent. Our agents are exogamous and only realise their potential for producing new agents when there is an empty cell within their fradius

As a default, the deathrate, or probability of an agent household disbanding in each time step, is set to the same value as the birthrate which means the net population growth rate is zero. The LGM is generally thought to represent a period of initial population contraction, followed by a period of little or no population growth in Europe, so assuming zero growth in our initial model is an appropriate first step [40]. We then experiment with variable effects of the environment on the deathrate as explained below.

We used two initial experiments (see SI) to ensure that our model produces patterns consistent with the growth rates reported for ethnographically or archaeologically described hunter-gatherer populations, which is between 0 and 0.04% (e.g., [43–46]). We then test the impact of habitat suitability on agent deathrates which in turn affects the size and distribution of regional populations in Europe. This has the dual advantage of allowing us to explain the components of the final model piece by piece and carefully investigate the effects of each additional component on the model results.

The Model: Impact of habitat suitability on survival

Experiments 1 and 2 held the growth rate at 0%, which meant setting the probability of disbanding (the “deathrate”) equal to the probability of creating a new household (the “birthrate”). In both experiments, HS-weighted random walks produced a spatial clustering of agents in high habitat suitability cells resulting in slightly larger populations in these areas. In the final model, we test the impact of habitat suitability (HS) on regional population structure and growth rates and so, in addition to weighting the random walks using HS, we scale the agents’ probability of survival to the HS value of their current cell at each time-step. We expect that groups of agents in more suitable habitats (i.e., regions with relatively high HS values) will experience higher growth rates due to decreased deathrates, while groups in more “risky”, marginal habitats (low HS values) will experience low to negative growth, with knock-on effects on regional population size and structure. There are no archaeological or ethnographic datasets to help us determine what the magnitude of the effect of HS on mortality and thus, population growth rates should be. In this experiment, therefore, we test different magnitudes of effect, observing their impact on overall growth rate. The null hypothesis is that the population’s growth rate should be within the limits described in the archaeological literature, that is between 0 and 0.04% [40,44].

In the final model we assume a simple, linear relationship between the probability of an agent persisting (deathrate) and HS. The HS index has a scale of 0 to 1, where 0 is highly unsuitable and 1 is highly suitable. We cannot simply assume that HS 1 = 100% probability of persisting and conversely, that HS 0 = 0% probability of persisting, however, because our null hypothesis is that the growth rate should approach 0% and the birthrate (probability of a new agent being created) is fixed to 17.6% / yr (see section 2.2.2, above). Under these conditions, a growth rate ~0% would only occur when HS values approach HS = 1 (which is extremely rare). As a result, the overall growth rate would be far lower than expected. A trial run confirmed that this setting results in total extinction of the metapopulation in a very short time. Therefore, we need to set up a new parameter that allows us to test the impact of HS on survival under conditions consistent with expected growth rates.

The new model parameter, hs_benefit, (Eq 1) creates the relationship: (1) where b is the birthrate, a is hs_benefit, d_h0.5 is the deathrate at HS 0.5 and d_h1.0 is the deathrate at HS = 1.0. We set the deathrate equal to the birthrate at HS = 0.5, which means zero growth occurs on cells at the middle of the HS value range. Using eq.1 we then derive the slope and y-intercept describing the linear relationship between HS and the deathrate allowing us to calculate the probability of persistence for all possible HS values (Fig 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. In experiments 1 & 2, habitat suitability had no effect on the deathrate which was set equal to the birthrate.

Here, we vary the benefit agents receive from being on higher HS values (and conversely, the penalty for being on cells with lower HS values) using the hs_benefit parameter.

https://doi.org/10.1371/journal.pone.0217996.g003

We use the hs_benefit parameter to test a range of possible relationships between the probability of persistence and habitat suitability in 25% increments. Concurrently, we conducted a series of runs with fradii values of 50, 100 and 200 km, using HS-weighted mobility. We conducted five runs of each parameter combination and each run was 3000 monthly time steps long, or 250 years. All experiments, including the code which exports the data tables and maps are included in the model code (https://doi.org/10.25937/na38-tj46). The processing time required for larger radii and greater hs_benefit runs is very high due to the large number of agents, as we shall see below. For these more computationally expensive runs, we shifted from our standard desktop processing to Graham, a high-performance computing system operated by the University of Waterloo.

The model results

At the beginning of the runs, the initial metapopulation of 7000 randomly distributed agents suffers a population decline as agents in areas distributed far from core areas suffer attrition. The metapopulation survives within the core regions, and then begins to grow from these core areas in most of the experiments (Fig 4). Note that all R scripts used to evaluate the model outputs and produce the figures below are available for download in our SocArXiv repository (http://doi.org/10.17605/OSF.IO/N24RQ).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Metapopulation size increase as a function of monthly time steps for five values of hs_benefit at three different scales of mobility.

https://doi.org/10.1371/journal.pone.0217996.g004

Fig 4 tracks the overall metapopulation size for different values of hs_benefit and of mobility (fradius) over the 3000 monthly time steps or 250 years of the runs. When hs_benefit = 0 the size of the population remains steady since HS does not affect the deathrate, which remains equal to the birthrate and so the growth rate is 0. As hs_benefit increases (Fig 4) the population grows at a rate roughly proportional to the magnitude of hs_benefit. In several of the experiments the metapopulation reaches the carrying capacity of the environment (~25000 agents) before the end of the run (Fig 4). In the other experiments (hs_benefit ≥ 0.25) the population was still growing when the run was terminated and would probably have reached the size threshold eventually. This reflects the fact that model parameters only allow agents to actualise the potential birthrate if they have available space for their “offspring” to move to, effectively treating HS as a proxy for carrying capacity. This has the effect of depressing the overall birthrate in regions with more spatially clustered resources (while maintaining the deathrate steady) independently of hs_benefit.

The model results show that if the structure of the environment (HS) has no impact on the probability of survival the initial population size remains stable. The distribution of the population is only weakly focussed on core areas. Increasing the chances of survival for agents at higher HS cell values (by increasing hs_benefit > 0) and conversely, decreasing the chance of survival on low HS cells, increases the size of the metapopulation as agents move onto cells with higher HS values (due to the weighted random walk) where their probability of survival is higher. The metapopulation grows until it reaches carrying capacity (Fig 4) though population growth is not uniformly distributed across the regions (Fig 5).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Final population size by “region” for different values of hs_benefit and fradius.

Region colors follow legend in the next figure and the map colors in Fig 2.

https://doi.org/10.1371/journal.pone.0217996.g005

The final population size of each region reflects the size of the core area, which contains higher HS values (Fig 3A) and, to a lesser extent, the size of the region (Fig 3B) and its average HS value. The Central (France) region, for example, has the largest core area, is the second largest region, and in addition to having the highest average HS score (Table 1), which results in a higher population growth relative to all other regions in all of the experiments. As regional populations grow, they expand from the core area towards the periphery of the region and will, where made possible by the radius of mobility, migrate into adjacent regions over generations of agent movements. The size of a region relative to its core area and the distance between core areas dictates the probability and direction of migration.

The genetic results (Fig 6) reflect the fact that regions with larger and better connected core areas have a greater chance of growth and subsequent migration. The well-connected Central region contains the largest population which, as it expands outward and migrates, usually comes to dominate the genetic profiles of other regional populations. Relatively higher rates of mortality, compared to the “birth” of new agents, outside the Central region due to lower average HS scores also favour the spread of “French” traits through migration and admixture. Nevertheless, the results also show that gene flow occurs between all of the regions, at all scales of mobility and under all test conditions except when mobility is constrained within a 50 km radius, and even then some gene flow occurs. At this lowest mobility condition of 50 km, Italy retains a larger proportion of its original traits reflecting its relatively isolated geographic position.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Genetic structure of regional populations for different values of hs_benefit and fradius.

https://doi.org/10.1371/journal.pone.0217996.g006

The weighted random-walk means that groups occupy much of the landscape occasionally, however, Figs 7 and 8 (below) show the frequency of occupation over the whole domain, which reflects the cumulative effect of higher population growth rates in the core regions. The Central region (the “France” region, which includes northeastern Spain) is the most intensively occupied region during the LGM and therefore reflects the location of a large “glacial refugium”. As higher hs_benefit values create a growing imbalance between low and high HS cells, regional disparities in occupation frequencies become more apparent, e.g. central Iberia and northern Europe are eventually abandoned in favour of the Central region when hs_benefit = 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7.

Relative occupation frequency for increasing degrees of hs_benefit with fradius = 50 km. hs_benefit values a) 0, b) 0.25, c) 0.50, d) 0.75, e) 1.0. n.b. Each map's color scheme is normalized to its own maximum occupation frequency value and is exported by the model to represent the final state of one run. European border polygons republished from GADM.org under a CC BY license, with permission from GADM, original copyright 2018.

https://doi.org/10.1371/journal.pone.0217996.g007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8.

Relative occupation frequency for increasing degrees of hs_benefit with fradius = 200 km. hs_benefit values a) 0, b) 0.25, c) 0.50, d) 0.75, e) 1.0.

https://doi.org/10.1371/journal.pone.0217996.g008

Implications of results

We reject the null hypothesis that habitat suitability does not affect rates of survival and thus population structure (hs_benefit = 0), because it produces demographic patterns at odds with the archaeological record. For example, it predicts a relatively large proportion of the population would inhabit the “Jura” region (Fig 5) which has a core area in the Jura but includes Central European countries as well as northern latitudes in Western Europe for which there is little or no evidence of human presence during the LGM (Fig 7).

If we accept the hypothesis that habitat suitability reflects a suite of environmental variables that affect human survival, however, we have to determine the magnitude of the effect. We tested different magnitudes of effect through the introduction of a parameter (hs_benefit) that modifies the probability of survival as a function of HS. In the absence of robust measures of productivity, it isn’t possible to calculate the carrying capacity of the different regions in our study but the habitat suitability (HS) index serves as a proxy. When HS affects the survival rate through incremental increases in hs_benefit, population growth increases until the carrying capacity of the environment is reached, at approximately 25,000 agents or 100,000 individuals, but remains within the upper range of archaeological estimates [33,40,44] in all of our experiments. Different magnitudes of hs_benefit are all possible, therefore, and the difference between them is mostly reflected in the speed with which the metapopulation reaches carrying capacity and the concentration of the population towards the core areas.

An hs_benefit value of 0.25 is the most “conservative” experiment and when human mobility is constrained to a radius of 50 km under these conditions, regional genetic signatures are evident and Italy retains a strong genetic identity. In subsequent experiments, with increasing penalties for occupying cells with low HS values, the metapopulation expands but the regional genetic signatures become more and more homogeneous at all levels of mobility, with the notable exception of Italy at the 50 km radius. At the same time, the level of human activity in regions beyond the core areas, e.g. northern Europe and central Iberia, diminishes (Fig 7 and Fig 8).

Under all of the test conditions where HS has an impact on human survival, the Central region (which includes France and northeastern Spain) has the largest population by far and several of the other regions struggle to maintain themselves. It is interesting to note that the Jura contains a fairly large population, which is not well reflected in the archaeological literature, although this region is receiving renewed attention [53–56]. Gene flow occurs between all of the regions, in all of the experiments, but genes associated with the Central “French” region dominate the regional genetic profiles. The dominance of the “French” genetic signature is due to its population size (which is determined by the size of its core area) and high connectivity with several adjacent core regions, which means that the expanding population has room to grow. This trend becomes more pronounced as the strength of the effect of HS on survivorship is increased, i.e., as the structure of the environment takes on a greater role in determining the chances of survival.

Our results suggest the important role the Central region would have played as a glacial refugium during the LGM. Its large core area is capable of sustaining a large population and leads to net population growth even as other regions experience more agent “deaths” than “births”. The strength of this effect increases as habitat suitability is given a greater role in determining the probability of survival (i.e. with increasing hs_benefit).

We interpret the results of these experiments in terms of a “source” and “sink” dynamic, similar to the one proposed for the Middle Pleistocene of Europe [57]. From this perspective, the Central region acts as the “source” of population growth, due to its centrality (connectedness), its large core area and its relatively high habitat suitability index. The other regions appear to have acted as “sinks”, with a net population loss being replenished by in-migration from the Central region. The genetic signature from “Valencia” (Mediterranean Spain), for example, disappears in two of our scenarios (Fig 6) albeit the least conservative of the hs_benefit scenarios. The source-sink dynamic has important implications for the genetic and cultural diversity of people living in Europe during the LGM.

We test three ranges of mobility (50 km, 100 km and 200 km). While lithic transfer patterns suggest mobility of up to 200 km, we suspect that movement on this scale would have been relatively rare, possibly falling into the category of “non-utilitarian” mobility described by Whallon [58]. Furthermore, with a radius of 200 km the implied territory size is twice as large as the highest value recorded for ethnographically known hunter-gatherer groups (Table 4–1 in [41]). It is possible that regions with more spatially dispersed patches of resources required higher levels of mobility, but the shorter ranges are probably better approximations of average mobility patterns during the LGM. A range of mobility of 50 km correspond to average lithic transfer distances recorded for Western Europe during the Upper Paleolithic [48] and produces annual territories of 7,854 km². This figure is compatible with the size of annual territories reported for hunter-gatherers who depend upon meat as an important resource [41]. Regional genetic variability is highest when mobility is restricted to a 50 km radius particularly in Italy, as one would expect given its poor connectivity with other regions. Since we did not include the Balkans in this study we cannot speak to possible gene flow from this region, however. When mobility is constrained to 50 km, land-use patterns concentrate around core areas and the “hinterlands” are abandoned, whereas at 100 and 200 km people still make use of these areas, although occupation intensity focusses on the core areas (compare Fig 8A and 8B). There is evidence in the archeological literature that regions outside of the traditional refugia, such as the Paris Basin [59,60] and Central Spain [61,62] could have been exploited during the LGM which could demonstrate the importance of at least occasional long-distance movements.

Genetics and material culture of the Last Glacial Maximum

Ancient DNA (aDNA) from human fossils dating to the LGM has the potential to validate or invalidate the conclusions we make in this paper. Indeed, our hypothesis (that HS affects the population dynamics of Western Europe during the LGM) and the ensuing models lead us to make very concrete predictions about regional population size, the genomic diversity and the structure of coalescent trees that researchers could construct from new samples taken from across European sites dating to the LGM. These include the existence of gene flow between all of the core areas within Western Europe, the dominance of genes from the Central region (including France and Northeastern Spain) which acts as a source population for Western Europe, and in some scenarios, the loss of genetic diversity from Mediterranean Spain (“Valencia”). In a scenario involving low human mobility patterns, Italy retains a distinct regional genetic signature.

Despite the increasingly rapid pace of sequencing, the LGM has been largely neglected so far [63] hampering our ability to test our predictions. Morphological data suggest that the human gene pool changed during the LGM [64]. Fu et al. [65] discuss genetic transitions surrounding the LGM at a continental scale, but most of their samples date to well before or well after the LGM and therefore cannot be used to evaluate the metapopulation structure during the LGM itself. Indeed, the only European genome of relevance is an individual from El Miron Cave, Spain dated to shortly after the LGM at ~18.7 ka and associated with the Magdalenian industry. Fu et al. [65] conclude, from comparing this genome to older European Gravettian genomes from ~34–26 ka and younger Magdalenian genomes from ~19–14 ka, that an endogenous European population survived the LGM within Europe and repopulated the continent afterwards. Similarly, Posth et al. [17] find a significant population contraction during the European LGM due to the loss of mtDNA haplogroup M in the post-LGM metapopulation, but they do not present genomes dating to the LGM itself. A greater number of ancient genomes, with finer spatial and temporal granularity, will be needed to test the predictions presented in this paper.

Our model does not include social or linguistic variables, nor do we examine cultural diffusion. However, our results have implications for cultural evolution. Proponents of the Cumulative Culture hypothesis have pointed out that demographic processes affect cultural transmission and rates of cultural innovation. For example, Powell et al. [4] show that levels of cultural innovation increase with the number and level of interaction between regional populations. Similarly, Derex and colleagues show that innovation rates may accelerate when the metapopulation is at an intermediate level of connectivity [66,67]. These studies and others suggest that the large and well-connected metapopulation we predict for the Central region will have favoured cultural transmission and made this region an important hub for cultural innovation. This prediction is compatible with the archaeological record. The Gravettian technocomplex which pre-existed the LGM was displaced by the Solutrean culture, which overlaps with the Early Magdalenian and Badegoulian in the Central region [68] and the Epigravettian in Italy [69]. The occupation frequencies predicted by the model do not suggest a direct connection between the core areas in the Iberian Peninsula (Atlantic and Mediterranean) but rather, suggest that these regions were indirectly connected through the central “French” region which may explain the cultural differences between these regions [18,70]. In the space of ~4,000 years, the Gravettian technocomplex, which was distributed across Western Europe prior to the LGM, is succeeded by the Solutrean in France and the Iberian Peninsula, which overlaps with the Early Magdalenian and the Badegoulian in France around 23,500 cal. BP [68]. The Epigravettian persists in Italy throughout the LGM [69].

Our results also predict that, unless we accept high levels of mobility and annual territories with radii of 100 km or 200 km, and high selective pressure (hs_benefit = 1), Italy forms a regionally distinct part of the Western European metapopulation, both in terms of its size and genetic composition. This observation is perfectly compatible with the material culture record, during the LGM the Epigravettian arose in Italy, distinguishing it from the other regions [69]. We will explore the dynamics of cultural interaction further in an upcoming paper.