The demographic response to Holocene climate change in the Sahara

The timing and development of Holocene human occupation in the now hyperarid Sahara has major implications for understanding links between climate change, demography and cultural adaptation. Here we use summed probability distributions from 3287 calibrated 14 C dates from 1011 archaeological sites to demonstrate a major and rapid demographic shift between 10,500 and 5500 years BP. This event corresponds with the African Humid Period (AHP) and is sub-continental in scale, indicating climate as the prime factor driving broad-scale population dynamics in northern Africa. Furthermore, by providing a high temporal resolution proxy for effective carrying capacity our population curve offers an inde- pendent estimate of environmental change in northern Africa, indicating a temporal delay in the terrestrial response to atmospheric climate change. These results highlight the degree to which human demography is a function of environment at the appropriate scale of observation in both time and space and sheds important new light on the social response to global environmental change. © 2014 The Authors. Published Elsevier


Introduction
The Sahara is the world's largest hot desert. With only 2.5 million people for over 3.5 million square miles, it currently supports one of the lowest population densities on earth, yet this has not always been the case. Over the course of the Holocene the Sahara underwent major climatic changes, marking the beginning of the African Humid Period (AHP) approximately 12,000 years ago (Gasse et al., 1990;Street-Perrott and Perrott, 1990;deMenocal et al., 2000;Adkins et al., 2006;McGee et al., 2013;Tierney and deMenocal, 2013), and what has been referred to as the 'Green Sahara' (Sereno et al., 2008;Drake et al., 2011;Dunne et al., 2012). Whilst the Holocene AHP is perhaps one of the most thoroughly documented and well-dated climate change events, its implications for ancient human settlement is less well understood.
Occupation is clearly testified in the frequent rock engravings that are scattered throughout the upland regions of the desert, illustrating a lush environment with Sahelian and riverine fauna and scenes of large-game hunting, livestock herding and religious ceremony (Mori, 1965;di Lernia and Gallinaro, 2010). Settlement studies from the Fezzan (Cremaschi and Zerboni, 2009;Cancellieri and di Lernia, 2013), Tenere (Sereno et al., 2008;Garcea, 2013) and Eastern Sahara (Wendorf et al., 2001;Hoelzmann et al., 2002;Kuper and Kr€ opelin, 2006), meanwhile, have shown evidence for intensive occupation, initially by pottery-using hunter-gatherers, and then by increasingly mobile pastoralists after the introduction of domestic livestock c. 7500 years BP (Marshall and Hildebrand, 2002;Dunne et al., 2012). Taken together, these lines of evidence strongly support a Holocene occupation of the desert, indicating close links between changing environments, human occupation and adaptation. However, although such regional perspectives provide vital information about the circumstantial influences of climate on Neolithic populations, little is known about the broad scale population trends, or the degree of synchronicity between demographic change and the onset and termination of the Holocene AHP.
Here we present results on the demographic reconstruction of the Saharan Holocene, revealing a temporal delay between the onset of humid conditions as represented in sedimentary dust flux records, and the human re-occupation of the desert. We estimate population growth and decline, and identify key stages in the Sahara's recent demographic history; a massive increase in trans-Saharan population levels shortly after 11,000 years BP, a temporary decline in population levels between 7600 and 6700 years BP, and a major population collapse at the end of the AHP, between 6300 and 5200 years BP.

Timing and abruptness of the Holocene African Humid Period
Debate is on going over the timing and abruptness of the AHP. Both sediment flux records from deep sea cores off the coast of north-west Africa (deMenocal et al., 2000;Adkins et al., 2006;McGee et al., 2013) and recent correlation analysis of dD wax isotopic values from east and northeast Africa (Tierney and deMenocal, 2013) point to a rapid hydroclimatic shift between 12,500e11,500 and 5500 years BP. Pollen and sedimentological records from Lake Yoa, in northern Chad however, indicate a more gradual decline in precipitation and subsequent progressive drying of the regional ecosystem (Kr€ opelin et al., 2008;Francus et al., 2013). This discrepancy is as yet unresolved, and is most likely a result of the differential sensitivity of the various proxies (Holelzmann et al., 2004;Francus et al., 2013;Tierney and deMenocal, 2013), although changes in regional hydroclimates may also have been amplified by vegetation feedback (Renssen et al., 2006), and local groundwater conditions (Holelzmann et al., 2004;L ezine et al., 2011;H ely et al., 2014). Patches of more arid terrain, for example, are hypothesized for parts of the Central Sahara based on pollen core records from the North Atlantic (Hooghiemstra et al., 1992) and Holocene lakes in the Central Sahara (H ely et al., 2014) as well as the fact that some sub-Saharan aquatic species do not exhibit trans-Saharan spatial distributions (Drake et al., 2011). A key issue therefore, is to assess the chronology of demographic change, and the spatial range over which those adaptations occurred.

Material and methods
In this paper we use a Monte-Carlo Summed Probability Distribution (MCSPD) method, which assesses the combined probability distribution of a large number of radiocarbon dates. We test the null hypothesis of an exponential demographic increase, which is based on the expectation of taphonomic loss in the archaeological record, and global exponential population increase across this broad timescale. Although there is increasing evidence to suggest the Sahara may have been more populated during the AHP, our null hypothesis must assume no prior knowledge of this in order to formally test the magnitude of demographic change, and whether the fluctuations in our SPD could be merely due to random sampling. The underlying principle of this method assumes a monotonic relationship between the population size and the amount of radiocarbon dates recovered, which is reliant on the law of large numbers to overcome small-scale temporal and spatial biases. Therefore by combining all 14 C dates in a region, the fluctuation in the density of dates through time can be used as a demographic proxy to investigate the timing of population change in northern Africa. Full details on the method can be found in Shennan et al. (2013) and Timpson et al. (in press), and all dates are calibrated using the IntCal09 curve (Reimer et al., 2009). Our initial analysis uses 3287 published radiocarbon dates from 1011 Neolithic archaeological sites (Fig. 1). Further details about the data and the full radiocarbon date list can be found in Appendix A.

Bootstrap analysis
We have also used a bootstrap analysis to provide an additional broad indication of the extent to which the Summed Probability Distribution (SPD) shape can be attributed to a genuine underlying pattern in the data, rather than random sampling fluctuations from a small sample. We generated 1000 bootstrapped SPDs by resampling (with replacement) the observed dates and calculated local 95% confidence intervals to indicate regions with greater confidence.

Calculating the effects of climate on human demography
Initially we attempted a regression analysis on our demographic curve and several high-resolution palaeoclimate proxies derived from sedimentary dust flux records (Adkins et al., 2006). However, as shown in Section 3.1 and 3.2, and discussed below, formal correlations were problematic. In order to test the hypothesis that climate change forced major population changes in the Saharan Holocene, we therefore divided our global dataset into four broad regions ( Fig. 1), the Atlas, Eastern Sahara, Central Sahara and the Western littoral and tested for correlation between the SPDs for these regions. The rationale behind this approach is that if populations respond concurrently across a diverse and large enough area, we can infer an exogenous cause such as climate to be the driving force. These sub-regions are purposefully very broad in order to maintain statistically viable sample sizes. Nonetheless, the boundaries between sub-regions are established in light of general geographic features, such as the distribution of palaeohydrological catchment areas, and to a lesser extent, cultural features such as the distinct pottery styles of the Western Littoral. Correlation coefficients between each of the four sub-regions are calculated using Pearson's R. However, a certain amount of similarity can be expected due to short-term fluctuations inherited from the calibration process, since all data is calibrated using the same calibration curve. Therefore p-values are estimated by comparing the correlation between the observed SPDs with the correlation between simulated SPDs. To ensure p-values are even more conservative, simulations are generated using a uniform null model, which fits the observed trends closer than exponential. Five hundred simulated SPDs are generated for each sub-region (as detailed in Shennan et al., 2013) and for each pair of sub-regions the correlation coefficient for all 250,000 permutations are calculated. The pvalue is then estimated as the proportion of these simulated Rvalues that are greater or equal to the observed R-value. All results are shown in Table 1.

Two-dimensional kernel density method for calculating population surface
Whilst SPD analysis performed on aggregated data provides high-resolution estimation of change through time, it provides no information on the spatial distribution of dates. Dividing the data into four geographically and culturally informed sub-regions provides four SPDs, which can be compared in order to draw broadscale inferences about differences between these regions. However, this approach has limitations since smaller sub-regions necessarily have smaller samples, so increasingly suffer from sampling Red represents the Eastern Sahara, blue the Central Sahara, orange the Atlas and purple the Western littoral. The grey dots represent sites that were not included in any of the four sub-regional analyses, but were included in the overall analysis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) error. Furthermore it becomes difficult to objectively justify where to draw the boundaries for increasingly arbitrary smaller subregions. Therefore, in order to more accurately analyse spatial variation we adopt a systematic approach, building on previous 2dimensional kernel density estimation methods (Collard et al., 2010), which does not require regional sub-division. Instead, the spatial and temporal location of each sample is considered as a point in space-time, and these points are then interpolated both through space and time. This is performed in discrete time-slices, resulting in many probability density surfaces, which are eventually stitched together to form a movie (Appendix B). Initially an SPD is generated for each of the 1011 archaeological sites (as detailed in the SPD method above). Smoothing through time is achieved by applying a rolling mean to each site's SPD so that for each time-slice (circa every 30 years) a much wider span of dates are included.
Smoothing through space is achieved using 2-dimensional kernel density estimation. One thousand and eleven kernels are generated (one for each site) using the spherical law of cosines to ensure kernels are perfect circles on the surface of a sphere. Each kernel is then weighted by its respective smoothed SPD. Smoothing parameters were selected heuristically, a 300 year rolling time bin for temporal smoothing, and a bandwidth of 3.7 km for spatial smoothing.

Results
Our analysis reveals a major demographic event with a significant departure from the null model (P < 0.0001) coinciding with the AHP, between 11,000 and 5500 years BP (Fig. 2). The first statistically significant period occurs around 10,250 years BP, albeit only for a brief time, followed by a second peak in population around 9500 years BP, and a more sustained growth trend until 7500 years BP. Between 7600 and 6700 years BP, the Sahara appears to have undergone a substantial population decline, which is consistent with records of occupational hiatus in the Tenere (Sereno et al., 2008). This is followed by a rapid recovery between 6700 and 6300 BP at which point population levels reach their relative optimum. Despite smaller scale fluctuations until 5500 years BP, this resurgence of Holocene populations appears to be only temporary, and is quickly followed by a major, and irreversed population collapse occurring between 6300 and 5200 years BP. Subsequent to this collapse, populations levels revert to a more gradual and permanent decrease albeit with a short reprieve between 3800 and 3500 years BP.

Climate as the prime factor driving broad-scale population dynamics in northern Africa
We were unable to detect any significant correlation between our SPD and the available palaeoclimate proxies. This may be due to variation in the chronological resolution of different palaeoclimate datasets, as well as their differing degrees of sensitivity (Francus et al., 2013;Tierney and deMenocal, 2013). However, it is also likely that these records do not accurately reflect the timing of AHP transitions across all of Africa. Whilst we have only included the sediment flux records from North West Africa it would certainly be worth incorporating additional datasets to test more regionally sensitive palaeoclimatic drivers. Our goal in this paper, however, is to examine the broad-scale relationship between long-term climate change and human demography. Therefore, in order to test the geographic extent of the population fluctuations, and hence the hypothesis that an exogenous force such as climate was the primary factor driving population dynamics, we divide the dataset into four sub-regions (Fig. 3). Analysis of the population curves shows a positive correlation for three of these four regions (Atlas, Eastern Sahara and Central Sahara), with R-values ranging from 0.62 to 0.92, and all significant at P < 0.0001 (Fig. 3, Table 1). These regions combine more than 3.4 million square miles, suggesting that climate was the driving force behind such a widespread demographic event. In particular, all three regions demonstrate synchronicity in the population increase at the start of the AHP, in the observed population fall and recovery between 7600 and 6300 years BP (both highlighted in red, Fig. 3), and more broadly (see Section 3.3) in the eventual decline around 6000 BP (not highlighted). In contrast, the Western littoral presents a very different demographic trend with no obvious concurrence with the AHP. The precise reason for this remains unclear, although the material culture suggests a distinct historical trajectory  . SPD for each of the four regional subdivisions. Three regions, which are significantly correlated, are drawn in black, whilst the Western Littoral, which does not correlate with any of the other regions, is drawn in grey. Thin lines represent the full SPD whilst thick lines show a 200-year rolling mean. The areas of particular synchronous interest are highlighted in red, i.e. the start of occupation around 11,000 years BP and the temporary population fall and rise between 7500 and 6500 years BP. The coloured dots on the map correspond to the sample key for each region as shown in Fig. 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (Vernet, 2007), with marine resource exploitation taking greater precedence over terrestrial resources, which may have provided greater resilience to climate forcing.

Delay in the terrestrial response to the onset of the AHP
By combining the 14 C dates from only the three broadly synchronous regions we were able to generate a refined SPD. Although this involved reducing the overall sample size, it filtered out the very different, and seemingly unrelated information from the Western Littoral (Fig. 4). This provides a more accurate and coherent population proxy to be utilized in more sensitive further analysis. Bootstrap analysis illustrates 95% confidence intervals in the refined SPD (grey ribbon Fig. 4), which is particularly narrow at the initial inferred demographic increase (shortly after 11000 years BP) and the subsequent decline (c. 6000e4500 BP).
Although the lack of formal correlation between our population curve and the available palaeoclimatic proxies may be partly due to differences in chronological resolution and differing degrees of sensitivity, suggesting that they may not encode all of the climatic variables that affect human populations. However, it may also be partly attributable to an apparent temporal delay (of approximately 1000 years) in the demographic response to the atmospheric changes characterizing the start of the AHP (Fig. 4). The first major population increase, which is evident in all three of our correlated sub-regions, occurs shortly after 11,000 BP, whilst the aeolian dust records indicate a shift in atmospheric conditions between 12,500 Table 1 Correlation coefficients and P-values for all four sub-regions.

Atlas
Central Sahara Eastern Sahara Central Sahara  and 11,500 BP (in fact deMenocal et al. (2000) propose an even earlier onset around 14,800 BP). To date, there is little consensus on the general terrestrial response to changing atmospheric conditions (Walther, 2010), although progressive response mechanisms on the scale of centuries to millennia have been hypothesized (Jones et al., 2009). Based on the spatio-temporal distribution of palaeohydrological features in the Sahara (L ezine et al., 2011), the first major phase of lake extension between 16 and 22 N occurs only around 10,500 years BP, and not until 9500 years BP north of 23 N, implying a delayed stabilization of the hydrological network some 1000 years after the increase in monsoon moisture flux inferred from the sediment records (Adkins et al., 2006). Similarly, our results indicate a corresponding delay in the demographic response, suggesting a progressive recovery of the terrestrial ecosystem, especially the plant and animal resources on which the early Holocene hunter-gatherers depended. In contrast, no such delay is observed at the population collapse and end of the AHP, with both the SPD and aeolian dust flux records, exhibiting remarkable synchronicity between 6300 and 5200 BP (blue polygon, Fig. 4).

Progressive population decline at the termination of the AHP
Following the onset of current arid conditions, between c. 6000 and 4500 BP, the demographic patterns across our three regions become less similar, suggesting divergent regional trajectories in the response to climatic deterioration. The Central Sahara, and to a lesser extent the Atlas region, appear to undergo more gradual decline in population density, as opposed to the abrupt population collapse observed in the Eastern Sahara. This may be partly explained by the emergence of pharaonic civilization along the Nile River leading to a spatial bias in our results (Wenke, 1991;Kuper and Kr€ opelin, 2006;Nicoll, 2012) as only Neolithic dates are included in this analysis. However, by this time, aridification was already considerably advanced in the eastern Sahara (Linst€ adter and Kr€ opelin, 2004;Tierney and deMenocal, 2013), and despite pharaonic occupation of the Nile corridor, the archaeological evidence does indicate abandonment of the desert regions (Kuper and Kr€ opelin, 2006). Runoff from the central Saharan massifs, meanwhile, would have played a major role in sustaining ground water levels in some central lowland regions (L ezine et al., 2011;Krinner et al., 2012;Francus et al., 2013), providing a temporary buffer for pastoral groups as new adaptive strategies (i.e. pastoral nomadism e see Section 4) emerged over the course of the later Holocene (di Lernia, 2006;Tafuri et al., 2006).

Estimating population growth and decline rates
Analysis of the refined SPD from the three synchronous regions further enabled us to estimate variation in growth and decline rates. The initial period 14,000e11,000 years BP (Fig. 4-a) shows a steady background increase averaging 1.1% per 25-year generation. Shortly after 11,000 years BP there is a large increase in population density with a massive increase in the growth rate to 7.7% per generation, resulting in an overall 143% increase in only 300 years between 10,900e10,600 BP ( Fig. 4-b), after which the growth rate slows to an average 0.8% increase for the next 3000 years (Fig. 4-c). Between 7600 and 6300 years BP we see a loss of 2.1% per generation resulting in an overall 34% loss in the 500 years from 7600 to 7100 BP, soon followed by rapid growth (Fig. 4-e) of 3.1% per generation, totalling a 63% increase in the 400 years from 6700 to 6300 BP. At this point population levels reach their relative maximum over the 12,000-year study period. The demographic collapse at the end of the AHP (Fig 4-f) resulted in a loss of 1.8% per generation, a 55% loss over 1100 years between 6300 and 5200 years BP, after which population levels revert to a more gradual and permanent decrease of 0.7% loss per generation.

Discussion
The refined SPD allows us to explore the temporal qualities of demographic change, however it provides no inherent information about spatial variation, and hence the context of regional demographic changes. This is less problematic at the broad scale of analysis, when investigating the sub-continental response to climate change, but creates a practical limit on the resolution at which regional heterogeneity can be investigated. We therefore apply a 2-dimensional kernel density estimation method (see Section 2.3) to generate a probability surface at discrete moments in time across the entire Saharan region.
On a spatial level, our results highlight a number of key stages in the demographic history of the Saharan Holocene (Fig. 5). Prior to 11,000 years BP (Fig. 5A), low-density population groups made initial incursions into the central Sahara, particularly in the Tibesti Mountains, where increased hydrological activity is observed from as early as 15,000 years BP (J€ akel, 1979; L ezine et al., 2011). The Er Rif region in northern Morocco is also well represented during this early phase with populations including both late Pleistocene Iberomaurusian groups and early Holocene hunter-gatherers. These groups can be subsumed under the term epipalaeolithic (Linst€ ader et al., 2012), and represent a distinctly North African trajectory, which appears to have peaked around 13,000 years BP with a marked intensification in cave use (Bouzouggar et al., 2008). After 11,000 years BP (Fig. 5B), increasing population levels are represented across much of the northern and eastern Sahara, reflecting the shift observed in our three regional analyses (Fig. 3). Also at this time, atmospheric conditions appear to have stabilized (Fig 4), allowing for the development of savannah-like environments into the previously hyper-arid Sahara. As animals began re-populating these regions (Drake et al., 2011), human groups appear to have rapidly followed, with several key regions being a focus for occupation i.e. the eastern Sahara and northern Air Mountains (Fig. 5C).
Over the next three millennia there is a progressive increase in the representation of early Holocene human groups across the Central and Eastern Sahara (Fig. 5CeE), suggesting a westward expansion. Despite a temporary demographic decline (Fig. 2) and extensification in the settlement patterns around 9000 years BP in the central Sahara, population levels reach a peak around 7500 BP (Fig. 5E). This intensification of key regions and expansion of populations into new Saharan territories corresponds with the development of the Khartoum Mesolithic and spread of the so-called 'Aqualithic' culture (Sutton, 1977). These pottery-using, huntergatherer-fisherfolk utilized a broadly similar toolkit consisting of microliths, harpoon points, bone fish hooks and comb impressed ceramics, all of which are present across much of the Sahara by 10,000 years BP. The technological and stylistic homogeneity of early Holocene material culture (Yellen, 1998) indicates a rapid process of population expansion, which was facilitated by the extensive network of lakes, rivers and inland deltas (Drake et al., 2011).
Following the temporary population collapse between 7500 and 6500 years BP, which is most apparent in the central Sahara (Fig. 5F), there is a renewed peak in population levels across all regions. Although domestic livestock are present in the eastern Sahara by at least 8000 years BP and tentatively in the Acacus by 7500 years BP, the widespread adoption of livestock pastoralism, did not occur until the mid 7th millennium BP (di Lernia, 2006). This correlates well with the subsequent boom we see in the demographic curve around 6250 years BP, and the expansion of populations into new ecological niches in the central Sahara Fig. 5. Still frames from the kernel density animation showing the key stages of demographic development across time and space. A shows initial Early Holocene incursions into the Sahara between 14,000 and 11,000 years BP, BeC show gradual population increase across the Sahara between 11,000 and 10,000 years BP, especially in the eastern Sahara and northern Air mountains, D shows exponential growth especially in the eastern Sahara, E shows the early Holocene population peak around 7500 years BP with most regions, except the western littoral being populated, F shows population retraction, especially in the central Sahara and the start of occupation along the western littoral, G shows renewed population expansion in the Middle Holocene, H shows retraction in all areas coinciding with the termination of the AHP, and I shows the de-population of the Sahara and southwestward movement of late Holocene populations. (Fig. 5G). Dairying practises are also in evidence at this time in the Acacus region (Dunne et al., 2012), alongside the development of a wider ideological phenomenon centred on domestic cattle and found across the eastern and central Sahara (Paris, 2000;di Lernia, 2006). It is unclear whether the rapid spread of domestic livestock around 7000 years BP was prompted by climatic events, although there is evidence for increasing aridity, and higher inter-annual variability in rainfall after 6500 BP (Kr€ opelin et al., 2008), which may explain the temporary population decline. However, in contrast to the early Holocene, the introduction of domestic livestock appears to have altered the level of resilience amongst Saharan populations. The development of pastoral nomadism as an adaptive mechanism for dealing with the onset of arid conditions meant that some populations in the central and southwestern Sahara persisted long after the termination of the AHP (Fig. 5I).

Conclusion
The analysis presented here offers a first insight into the demographic response to Holocene climate change in the Sahara. Over a relatively short period of time climatic amelioration allowed human populations to colonize a vast region of the African continent, and yet within less than 6000 years, they were gone. More generally, this analysis offers greatly improved temporal resolution on effective carrying capacity, providing critical insight into the terrestrial response to changing atmospheric conditions, and a complementary proxy for the timing of Holocene climate change in northern Africa. We find evidence for a temporal delay of around 1000 years in the terrestrial recovery of the Saharan ecosystem at the start of the Holocene AHP. Our results also highlight a substantial and broad-scale demographic decline between 7500 and 6500 years BP, which is in accordance with major cultural and economic changes, but has yet to be identified in any current palaeoclimate proxies although recent work by H ely et al. (2014) suggests a possible synchronous dip in Saharan biodiversity, providing a potential source of corroboration. Further work is needed to explore the implications of this event, nevertheless it seems likely, based on the broad scale synchronicity between our sub-regions, that an exogenous force such as climate must have played a role.