Homo sapiens in Arabia by 85,000 years ago.

Understanding the timing and character of the expansion of Homo sapiens out of Africa is critical for inferring the colonization and admixture processes that underpin global population history. It has been argued that dispersal out of Africa had an early phase, particularly ~130-90 thousand years ago (ka), that reached only the East Mediterranean Levant, and a later phase, ~60-50 ka, that extended across the diverse environments of Eurasia to Sahul. However, recent findings from East Asia and Sahul challenge this model. Here we show that H. sapiens was in the Arabian Peninsula before 85 ka. We describe the Al Wusta-1 (AW-1) intermediate phalanx from the site of Al Wusta in the Nefud desert, Saudi Arabia. AW-1 is the oldest directly dated fossil of our species outside Africa and the Levant. The palaeoenvironmental context of Al Wusta demonstrates that H. sapiens using Middle Palaeolithic stone tools dispersed into Arabia during a phase of increased precipitation driven by orbital forcing, in association with a primarily African fauna. A Bayesian model incorporating independent chronometric age estimates indicates a chronology for Al Wusta of ~95-86 ka, which we correlate with a humid episode in the later part of Marine Isotope Stage 5 known from various regional records. Al Wusta shows that early dispersals were more spatially and temporally extensive than previously thought. Early H. sapiens dispersals out of Africa were not limited to winter rainfall-fed Levantine Mediterranean woodlands immediately adjacent to Africa, but extended deep into the semi-arid grasslands of Arabia, facilitated by periods of enhanced monsoonal rainfall.

Background

Homo sapiens evolved in Africa in the late Middle Pleistocene1. Early dispersals out of Africa are evidenced at the Levantine site of Misliya at ~194-177 ka2, followed by Skhul and Qafzeh, where H. sapiens fossils have been dated to ~130-100 and ~100-90 ka respectively3. While the Levantine fossil evidence has been viewed as the onset of a much broader dispersal into Asia4–6, it has generally been seen as representing short-lived incursions into the woodlands of the Levant immediately adjacent to Africa, where relatively high precipitation is produced by winter storms tracking across the Mediterranean7,8. While the Levantine record indicates the subsequent local replacement of early H. sapiens by Neanderthals, the failure of early dispersals to extend beyond the Levant is largely inferred from interpretations of genetic data9. Genetic studies have suggested that recent non-African populations stem largely10, if not entirely9, from an expansion ~60-50 ka, but this model remains debated. The absence of low latitude Pleistocene human DNA and uncertainties regarding ancient population structure undermine conclusions drawn from genetic studies alone. The paucity of securely dated archaeological, palaeontological and ancient DNA data - particularly across southern Asia - has made testing dispersal hypotheses challenging4,7,11.

Recent fossil discoveries in East Asia indicate that the early (particularly Marine Isotope Stage 5) dispersals of Homo sapiens extended across much of southern Asia. At Tam Pa Ling in Laos, Homo sapiens fossils date to between 70 and 46 ka12. Teeth assigned to Homo sapiens from Lida Ajer cave, Sumatra, were recovered from a breccia dating to 68 ± 5 ka, with fauna from the site dating to 75 ± 5 ka13. Several sites in China have produced fossil material claimed to represent early Homo sapiens14. These include teeth from Fuyan Cave argued to be older than 80 ka based on the dating of an overlying speleothem a few metres from the fossils15, and teeth from Luna Cave that were found in a layer dating to between 129.9 ± 1.5 ka and 70.2 ± 1.4 ka16. Teeth and a mandible from Zhiren Cave, China, date to at least 100 ka and have been argued to represent Homo sapiens, but other species attributions are possible17. The recent documentation of a human presence in Australia from ~65 ka is consistent with these findings18. Likewise, some interpretations of genetic data are consistent with an early spread of Homo sapiens across southern Asia10. These discoveries are leading to a radical revision of our understanding of the dispersal of Homo sapiens, yet there remain stratigraphic and taxonomic uncertainties for many of the east Asian fossils14,19, and thousands of kilometers separate these findings from Africa.

The Arabian Peninsula is a vast landmass at the crossroads of Africa and Eurasia. Growing archaeological evidence demonstrates repeated hominin occupations of Arabia20,21 each associated with a strengthened summer monsoon which led to the re-activation of lakes and rivers22–24, as it did in North Africa25. Here we report the discovery of the first pre-Holocene human fossil in Arabia, Al Wusta-1 (AW-1), as well as the age, stratigraphy, vertebrate fossils and stone tools at the Al Wusta site (Fig. 1, see also Supplementary Information).

Figure 1

Al Wusta location, map of site and stratigraphy.

A: The location of Al Wusta and other key MIS 5 sites in the region11; B: Al Wusta digital elevation model showing location of AW-1 phalanx, marl beds, lithics and vertebrate fossils, and the locations of the trenches and sections. The inset shows a satellite image of the site; C: Stratigraphic log of Al Wusta showing the sedimentology of the exposed carbonate beds, isotopic values, OSL ages for sand beds and U-series and ESR ages for AW-1 and WU-1601. Sands are shown in yellow: lower massive sands are aeolian (Unit 1), upper laminated sands are waterlain (Unit 3a) and have been locally winnowed to generate a coarse desert pavement (Unit 3b), lacustrine marls are shown (Unit 2) in grey (for full key and description see Supplementary Figures 13 and 14 and Supplementary Information 5). Section PD40 is shown as it contains the thickest sequence and is most representative of Al Wusta, chronometric age estimates (marked *) from the site are depicted in their relative stratigraphic position, see Supplementary Figure 14 for their absolute positions.

Results

AW-1 is an intermediate manual phalanx, most likely from the 3rd ray (Fig. 2a, Supplementary Information 1: see below for detail on siding and species identification). It is generally well-preserved, although there is some erosion of the cortical/subchondral bone, and minor pathological bone formation (likely an enthesophyte) affecting part of the diaphysis (Supplementary Information 1). The phalanx measures 32.3 mm in proximo-distal length, and 8.7 mm and 8.5mm in radio-ulnar breadth of the proximal base and midshaft, respectively (Supplementary Table 1).

Figure 2

Photographs and micro-CT scans of Al Wusta-1 Homo sapiens phalanx.

A: photographs in (left column, top to bottom) distal, palmar and proximal views, and (middle row, left to right) lateral 1, dorsal and lateral 2 views. Micro-CT cross-sections (illustrated at 2x magnification) include B (54% from proximal end) and C (illustrating abnormal bone).

AW-1 is more gracile than the robust intermediate phalanges of Neanderthals26–28, which are broader radio-ulnarly relative to their length and have a more 'flared' base. AW-1’s proximal radio-ulnar maximum breadth is 14.98 mm, which provides an intermediate phalanx breadth-length index (proximal radio-ulnar maximum breadth relative to articular length) of 49.6. This is very similar to the mean (± SD) for the Skhul and Qafzeh H. sapiens of 49.7 (± 4.1) and 49.1 (± 4.0) for Upper Palaeolithic Europeans, but 1.89 standard deviations below the Neanderthal mean of 58.3 (± 4.6)29.

To provide a broad interpretive context for the Al Wusta phalanx, we conducted linear and geometric morphometric (GMM) landmark analyses (Supplementary Information 1) on phalanges from non-human primates, fossil hominins and geographically widespread recent H. sapiens. Comparative linear analyses (Supplementary Information 1, Supplementary Tables 2 and 3, Supplementary Figure 1) reveal that there is substantial overlap across most taxa for all shape ratios, so AW-1 falls within the range of variation of H. sapiens, cercopiths, Gorilla, Australipithecus afarensis, A. sediba and Neanderthals. However, AW-1 is most similar to the median value or falls within the range of variation of recent and early H. sapiens for all shape ratios.

Geometric morphometric (GMM) analyses of AW-1 and various primate groups including hominins (see Supplementary Table 4 and Supplementary Figure 2 for landmarks, and Supplementary Table 5 for sample) are illustrated in Figure 3 and Supplementary Figure 3. PC1 and PC2 together account for 61% of group variance in shape. AW-1 is separated on these two shape vectors from the non-human primates and most of the Neanderthals. AW-1 falls closest to the recent and early H. sapiens and is clearly differentiated from all non-human primates. This is also shown by the Procrustes distances from AW-1 to the mean shapes of each taxonomic group (Supplementary Table 6).

Figure 3

Scatterplot of the first two principal components (PC) scores of the geometric morphometric analysis of the Al Wusta-1 phalanx compared with a sample of primates, including hominins.

Non-human hominoids: lilac; Gorilla: circles, Pan: triangles. Cercopithecoids: red; Colobus: triangles, Mandrillus: squares, Papio: circles. Neanderthals: blue diamonds. H. sapiens: green; early H. sapiens: circles, Holocene H. sapiens: squares. Al Wusta-1: black star, circled in red.

Three of the Neanderthal phalanges (from Kebara 2 and Tabun C1) are quite disparate from the main Neanderthal cluster and fall closer to the H. sapiens and Al Wusta cluster on PC1 and 2 (Figure 3 and Supplementary Figure 3). Having established the hominin affinity of AW-1, shape was analysed in more detail using a smaller hominin sample for which ray number and side were known, which included Kebara 2 and Tabun C1. The broader primate sample used in the first GMM analysis was not used for the more detailed shape analysis, as the initial comparisons show clearly that AW-1 is not a non-human primate and including this level of variation could potentially mask more subtle shape differences between hominins. The side and ray are also not known for most of the Neanderthal and non-human primate samples, meaning it would be impossible to evaluate the effect of these factors using this sample.

The more in-depth shape comparison and modelling using the hominin sample of phalanges of known ray and side (Supplementary Table 7) demonstrates that the long and slender morphology of AW-1 falls just outside the range of variation of comparative Middle Palaeolithic modern humans, but that its affinity is clearly with H. sapiens rather than Neanderthals (Fig. 4, Supplementary Table 8). Although both Pleistocene H. sapiens and Neanderthal landmark configurations fall almost completely inside the scatter for the Holocene H. sapiens sample in the principal components analysis (Figure 4), AW-1 is closest to Holocene H. sapiens 3rd intermediate phalanges. AW-1 overlaps with the Holocene H. sapiens sample, but is separated from the Pleistocene H. sapiens specimens by a higher score on PC2 and from the Neanderthal group by a simultaneously higher score on PC1 and PC2. The Procrustes distances (Supplementary Table 8), also show that AW-1 is most distinct from the Neanderthal phalanges, which fall towards the lower ends of both PCs and are characterised by shorter and broader dimensions. PC1 and PC2 in this analysis show that AW-1 is taller and narrower (in all directions: dorso-palmarly, proximo-distally and radio-ulnarly) than almost all the phalanges in the comparative sample and is particularly distinct from most of the Neanderthal phalanges. In this analysis AW-1 is closest in shape to 3rd phalanges of individuals from (in descending order of proximity) Egyptian Nubia, and Medieval Canterbury (UK), and Maiden Castle (Iron Age Dorset, UK) (Supplementary Table 9), although there is not a great difference in its distance to any of these specimens. These analyses suggest that the AW-1 phalanx is likely to be a 3rd intermediate phalanx from a H. sapiens individual.

Figure 4

Scatterplot of the first two principal component (PC) scores from the geometric morphometric analyses of AW-1 and sample of comparative hominin 2nd, 3rd, and 4th intermediate phalanges.

Wireframes show mean configuration warped to extremes of PC axes in dorsal (left), proximal (middle) and lateral (right) views. Convex hulls added post-hoc to aid visualisation.

The third ray is the most symmetrical ray in the hand and is therefore difficult to side, particularly when not all of the phalanges of a particular individual are present. Comparing AW-1 separately to right and to left phalanges (Supplementary Information 1.4) gives results which are very similar to the pooled sample, such that AW-1 is closest to Holocene H. sapiens 3rd rays for both right and left hands (Supplementary Figure 4, Supplementary Table 10). There is little difference in morphological closeness between AW-1 and its nearest neighbour in the samples of right and left bones (Supplementary Table 11), reflecting the lack of difference in morphology between the sides. It is therefore not possible to suggest whether AW-1 comes from a right or a left hand using these analyses.

AW-1 is unusual in its more circular midshaft cross-sectional shape (Fig. 2B), which is confirmed by cross-sectional geometric analyses (Supplementary Information 1.5). This may reflect the pronounced palmar median bar that makes the palmar surface slightly convex at the midshaft rather than flat, the latter being typical of most later Homo intermediate phalanges. However, more circular shafts may reflect greater loading of the bone in multiple directions and enthesophytes are a common response to stress from high levels of physical activity30. This morphology may reflect high and varied loading of the fingers during intense manual activity.

To determine the age of AW-1, and associated sediments and fossils, we used a combination of uranium series (U-series), electron spin resonance (ESR) and optically stimulated luminescence (OSL) dating (Methods, Supplementary Information 2 and 3). U-series ages were produced for AW-1 itself (87.6 ± 2.5 ka) and hippopotamus dental tissues (WU1601), which yielded ages of 83.5 ± 8.1 ka (enamel) and 65.0 ± 2.1 ka (dentine). They should be regarded as minimum estimates for the age of the fossils. In addition, a combined U-series-ESR age calculation for WU1601 yielded an age of 103 +10/-9 ka. AW-1 was found on an exposure of Unit 3b, and WU1601 excavated from Unit 3a, one metre away (Fig 1b). Unit 1 yielded OSL ages of 85.3 ± 5.6 ka (PD17), 92.2 ± 6.8 ka (PD41) and 92.0 ± 6.3 ka (PD15), while Unit 3a yielded an OSL age of 98.6 ± 7.0 ka (PD40). The OSL age estimates agree within error with the US-ESR age obtained for WU-1601 and the minimum age of ~88 ka obtained for AW-1. These data were incorporated into a Bayesian sequential phase model31 which indicates that deposition of Unit 1 ceased 93.1 ± 2.6 ka (Phase 1: PD15, 17, 41) and that Units 2 and 3 and all associated fossils were deposited between 92.2 ± 2.6 ka and 90.4 ± 3.9 ka (Phase 2: all other ages) (Supplementary Information 4, Supplementary Figure 11).

This ~95-86 ka timeframe is slightly earlier than most other records of increased humidity in the region in late MIS 5 32,33, which correlate with a strengthened summer monsoon associated with an insolation peak at 84 ka (Fig. 6). The underlying (Unit 3) aeolian sand layer at Al Wusta correlates with an insolation minimum at the end of MIS 5c. The chronometric age estimates for the site suggest that lake formation and the associated fauna and human occupation occurred shortly after this in time. Regional indications of increased humidity around the 84 ka insolation peak include speleothem formation at ~88 ka in the Negev34, and the formation of sapropel S3 beginning ~86 ka35. In both the Levant and Arabia, records are consistent with this switch from aridity to humidity around this time32–40. Precisely reconstructing regional palaeoclimate at this time and relating it to human demographic and behavioural change has proved challenging. This reflects both rapid changes in climate, as well as the complexities involved in dating relevant deposits41. In summary, combining chronological data (Supplementary sections 2-4), interpretation of the sedimentary sequence (described below), and the regional setting of Al Wusta, we conclude that lake formation and associated finds such as the AW-1 phalanx relate to the late MIS 5 humid period associated with the 84 ka insolation peak.

Figure 6

The chronological and climatic context of Al Wusta.

The Al Wusta lake phase falls chronologically at the end of the time-range of MIS 5 sites from the Mediterranean woodland of the Levant (~130-90 ka) and earlier than the late dispersal(s) (~60-50 ka) as posited in particular by genetic studies. The chronology of these dispersals and occupations correspond with periods of orbitally modulated humid phases in the eastern Mediterranean36 that are important intervals for human dispersals into Eurasia, and are also proposed to correspond with episodes of monsoon driven humidity in the Negev and Arabian desert34. Environmental amelioration of the Saharo-Arabian belt, therefore, appears to be crucial for allowing occupation at key sites that document dispersal out of Africa. A: East Mediterranean speleothem δO record from Soreq and Pequin Caves36; B: global δO record37; C: Insolation at 30 degrees north38, showing the temporal position of key sites relating to dispersal out of Africa2,3,11,48. The chronology for Al Wusta shows the phases defined by the Bayesian model at 2σ.

The sedimentary sequence at Al Wusta consists of a basin-like deposit of exposed carbonate-rich sediments (Unit 2, 0.4-0.8 m thick), underlain by wind-blown sand (Unit 1) and overlain by water-lain sands (Unit 3). The carbonate rich sediments of Unit 2 are interpreted as lacustrine marl deposits on the basis of their sedimentology, geochemistry, and diatom palaeoecology (Figure 1c, Methods, Supplementary Information 5). At both the macro- and micro-scale, these beds are relatively massive and comprise fine-grained calcite, typical of material precipitating and accumulating in a still-water lacustrine environment42. At the micro-scale there is no evidence for the desiccation or fluctuation of water levels typical of palustrine/wetland environments42, implying that the lake body was perennial. The diatom flora support this, containing species such as Aulacoseira italica and Aulacoseira granulata throughout the sequences, indicating an alkaline lake a few metres deep. The water was fresh, not saline or brackish, since saline tolerant species and evaporitic minerals are absent throughout. While δ18O and δ13C values of continental carbonates are controlled by a wide-range of variables, the values derived from the Al Wusta marl beds are compatible with the suggestion of marls precipitated in a perennial lake basin. The Al Wusta carbonate beds therefore indicate a perennial lake body a few metres in depth. The existence of a marl precipitating lake basin implies that this system was groundwater fed (to allow for sufficient dissolved mineral material to be present in the lake waters). Although the Al Wusta sequence represents a single lake basin, the development of such a feature over highly permeable aeolian sands in a region where no lake systems exist at the present day implies a local increase in water table that would require an increase in mean annual rainfall. Consequently, the Al Wusta sequence represents the occurrence of a humid interval at this time. The Unit 2 marl is overlain by a medium-coarse sand (Unit 3) with crude horizontal laminations, occasional clasts, fragments of ripped up marl and shells of Melanoides tuberculata and Planorbis sp. While some vertebrate fossils and lithics were found in the upper part of Unit 2, most were found in or on the surface of Unit 3. Unit 3a sands are waterlain and represent the encroachment of fluvial sediment as the lake environment shallowed and contracted. Unit 3b represents a winnowed lag formed by aeolian deflation of 3a. The sequence is capped by a dense network of calcitic rhizoliths marking the onset of fully terrestrial conditions.

A total of 860 vertebrate fossils were excavated from Unit 3 and the top of Unit 2 (n=371) and systematically surface collected (n=489). These include specimens attributed to Reptilia, Aves, and Mammalia (Supplementary Table 19, Methods, Supplementary Information 6). Notable taxa now extinct in Arabia are predominately grazers and include Hippopotamus, Pelorovis, and Kobus. The faunal community demonstrates a clear preference for temperate to semi-arid grasslands, and the presence of Hippopotamus and Kobus indicate permanent muddy, fluvial, or lacustrine conditions43 not currently found in the Nefud Desert, but consistent with the geological evidence from the site. The faunal assemblages show a strong affinity to African fauna, particularly Hippopotamus, Pelorovis, and Kobus44. Many large tooth pits on fossils indicate that large carnivores played a role in the accumulation of the deposit. Long bone circumference, completeness and numbers of green fractures suggests modification of bones by bone-breaking agents such as large carnivores or hominins (Supplementary Information 6). However, no evidence of cut-marks or hammerstone damage to the bones was observed.

An assemblage of 380 lithic artefacts (stone tools) was recovered from the excavation of upper Unit 2 and Unit 3 and systematic surface collection (Methods, Figure 5, Supplementary Information 7). They are of Middle Palaeolithic character and most are chert and quartzite. The assemblage demonstrates a focus on centripetal Levallois reduction, and is similar to other late Marine Isotope Stage 5 assemblages in the west and north of Arabia45, and contemporaneous assemblages in east (e.g. Aduma, BNS at Omo Kibish) and northeast Africa (e.g. Bir Tarfawi), as well as those from the Levant (e.g. Qafzeh)11 (Fig. 5).

Figure 5

Selected Al Wusta lithic artefacts.

A: argillaceous quartzite flake; B: quartz hammerstone; C: ferruginous quartzite Levallois flake; D: chert Levallois flake; E: Quartz recurrent centripetal Levallois core; F: quartzite preferential Levallois core with centripetal preparation and pointed preferential removal.

Discussion

Al Wusta-1 is the oldest directly dated H. sapiens fossil outside Africa and the Levant. It joins a small but growing corpus of evidence that the early dispersal of H. sapiens into Eurasia was much more widespread than previously thought. The site of Al Wusta is located in the Nefud desert more than 650 km southeast of Skhul and Qafzeh (Fig. 1A). This site establishes that H. sapiens were in Arabia in late MIS 5, rather than being restricted to Africa and the Levant as suggested by traditional models (Fig. 6). With Skhul dating to ~130-100 ka, Qafzeh to ~100-90 ka3,46 and Al Wusta to ~95-85 ka it is currently unclear if the southwest Asian record reflects multiple early dispersals out of Africa or a long occupation during MIS 5. The association of the Al Wusta site with a late MIS 5 humid phase (Fig. 6), suggests that significant aspects of this dispersal process were facilitated by enhanced monsoonal rainfall. While changes in behaviour and demography are crucial to understanding the dispersal process, climatic windows of opportunity were also key in allowing H. sapiens to cross the Saharo-Arabian arid belt, which often constituted a formidable barrier24,25.

Conclusion

Al Wusta shows that the early, Marine Isotope Stage 5, dispersals of H. sapiens out of Africa were not limited to the Levantine woodlands sustained by winter rainfall, but extended deep into the Arabian interior where enhanced summer rainfall created semi-arid grasslands containing abundant fauna and perennial lakes. After long being isolated in Africa1,47,48, the Late Pleistocene saw the expansion of our species out of Africa and into the diverse ecologies of Eurasia. Within a few thousand years of spreading into Eurasia our species was occupying rainforest environments and making long sea crossings to remote islands13,18. Adapting to the semi-arid conditions of the Saharo-Arabian arid belt represented a crucial step on this pathway to global success and the Al Wusta Homo sapiens fossil demonstrates this early ability to occupy diverse ecologies which led to us becoming a cosmopolitan species.

Methods

Site identification, survey and excavation

The site of Al Wusta (field code WNEF16_30) was discovered in 2014 as part of a programme of joint survey fieldwork of the Palaeodeserts Project, the Saudi Commission for Tourism and National Heritage, and the Saudi Geological Survey. It is located in the western Nefud desert, a few kilometres from the Middle Pleistocene fossil locality of Ti’s al Ghaddah49. The locations of all materials of interest (fossils, stone tools, geomorphological features, excavations and sample points) were recorded using a high-precision Trimble XRS Pro Differential GPS system and a total station, and entered into a GIS (Fig. 1). Elevation data (masl) were recorded as a series of transects across the site, and a digital elevation model (DEM) and contours interpolated (Spline) from all data with precisions of better than 10 cm in all (x,y,z) dimensions (22,047 points). This allowed visualisation and recording of the spatial relationships between materials in three dimensions (Fig. 1). Eight trenches were excavated into the fossil and artefact bearing deposits. These trenches revealed vertebrate remains and lithics, but no further human fossils were recovered.

Morphological analysis of Al Wusta-1 phalanx

The phalanx was scanned using micro-computed tomography (micro-CT) on the Nikon Metrology XT H 225 ST High Resolution scanner and X-Tek software (Nikon Metrology, Tring, UK) housed in the Cambridge Biotomography Centre, University of Cambridge, UK. Scan parameters were: a tungsten target; 0.5 mm copper filter; 150 kV; 210 mA; 1080 projections with 1000 ms exposure, and resulted in a voxel size of 0.02 mm3. The micro-CT data were reconstructed using CT-PRO 3D software (Nikon Metrology) and exported as an image (.tif) stack. Other CT data were obtained from the institutions cited in Supplementary Table 5 with permissions following the memoranda of understanding with each institution.

3D landmarks and semilandmarks were chosen to best describe the overall shape of the morphology of the AW-1 phalanx (Supplementary Table 4, Supplementary Figure 2), and were digitised on virtual reconstructions of phalanges created from micro-CT data in AVIZO 8 and 9.1 (FEI Software, Burlington, Mass.). Landmark coordinates were exported for use in Morphologika50. In Morphologika, generalized Procrustes analyses were performed to superimpose landmark coordinate data, and principal components analyses (PCA) were run to investigate similarities in shape between specimens. Shape differences along principal componentss were visualised and wireframes were produced in Morphologika, PC scores were exported to create graphs in R51. Procrustes distances between specimens were calculated using MorphoJ52. To avoid representing the same phalanges from different sides of a single individual as independent data points and to maximise sample sizes in pooled analyses, right phalanges were used in cases where the phalanges from both sides were present. Where only the left was present, this was used and ‘reflected’ (i.e. mirrored) in Morphologika to generate landmark configurations consistent with right phalanges.

U-series and combined US-ESR dating of fossil bone and teeth

The AW-1 phalanx (lab number 3675) and a hippopotamus tooth fragment (lab number WU1601) were collected from Trench 1 (Fig.1) for U-series and combined US-ESR dating, respectively. The external dose rate utilised the data of OSL sample PD40, which was collected in an equivalent position within unit 3a.

U-series analysis

U-series analyses were conducted at the Research School of Earth Sciences, The Australian National University, Canberra. The experimental setup for the U-series analysis of the phalanx was described in detail by Grün and colleagues53 (Supplementary Figures 2 and 3, Supplementary Information 2). Laser ablation (LA) was used to drill a number of holes into AW-1 following the approach of Benson and colleagues54. After a cleaning run with the laser set at a diameter of 460 μm, seven holes were drilled for 1000 s with the laser set at 330 μm. The isotopic data streams were converted into 230Th/234U and 234U/238U activity ratios and apparent Th/U age estimates and subsequently binned into 30 successive sections (each containing 33 cycles) for the calculation of average isotopic ratios and ages. A similar experimental setup and methodology were employed for the LA U-series analysis of tooth sample WU1601. The whole closed system U-series analytical datasets of the enamel and dentine sections were integrated to provide the data input for the ESR age calculations.

Combined US-ESR dating of the fossil tooth: ESR dose evaluation

The ESR dose evaluation of the hippo tooth was carried out at CENIEH, Burgos, Spain, following a similar procedure to that described in Stimpson and colleagues49. Enamel was collected from WU1601 and powdered <200 µm. The sample was then divided into 11 aliquots and gamma irradiated with a Gammacell-1000 Cs-137 source to increasing doses until 3.4 kGy. ESR measurements were carried out at room temperature with an EMXmicro 6/1 Bruker ESR spectrometer coupled to a standard rectangular ER 4102ST cavity. ESR intensities were extracted from T1-B2 peak-to-peak amplitudes of the ESR signal of enamel. Fitting procedures were carried out with a single saturating exponential (SSE) function through the pooled ESR experimental data derived from the repeated measurements, with data weighting by the inverse of the squared ESR intensity (1/I2) and following the recommendations by Duval and Grün55. Full details about the experimental conditions and analytical procedure may be found in Supplementary Information 2.

Combined US-ESR dating of the fossil tooth: Dose rate evaluation and age calculations

The combined US-ESR age of WU1601 was calculated with the DATA programme56 using the US model defined by Grün and colleagues57. The following parameters were used for the dose rate evaluation: an alpha efficiency of 0.13 ± 0.0258, Monte-Carlo beta attenuation factors from Marsh59, dose-rate conversion factors from Guerin and colleagues60, external sediment (beta and gamma) dose rate from the OSL sample PD40, a depth of 25 ± 10 cm, resulting in an age of 103 + 10/-9 ka.

Optically Stimulated Luminescence Dating

Three samples (PD15, PD17 and PD41) were collected from the aeolian sands (Unit 1) underlying the southern marl outcrop (Unit 2, Fig 1B). A fourth sample (PD40) was taken from the main fossil bearing bed (Unit 3). Individual quartz grains were measured on a Risø TL/OSL-DA-15 instrument using the single-aliquot regenerative-dose (SAR) method61. The burial dose for each sample (Db) was calculated using the central age model (CAM)62.

Environmental dose rates were determined using a Risø GM-25-5 low-level beta counting system63 (beta dose rate), field gamma spectrometry (gamma dose rate), and an estimate of the cosmic dose rate derived using site location and present day sediment burial depths64. Full optically stimulated luminescence dating methods and results are presented in Supplementary Information Section 3. All analyses were carried out in the Royal Holloway Luminescence Laboratory by SA and R C-W.

Age modelling

Chronometric ages for samples from the Al Wusta site were incorporated into a Bayesian sequential phase model implemented in OxCal v4.231 (Supplementary Information 4; Supplementary Figure 11. The model consists of two discrete phases separated by a hiatus. Phase 1 was defined by the three OSL ages (PD15, 17 and 41) for samples from the aeolian sands (Unit 1) underlying the lacustrine marls (Unit 2). Phase 2 was defined by the ages for the sand (PD40) and fossils (AW-1 and WU1601) from the waterlain sediments (Unit 3) overlying Unit 2. U-series ages for WU1601 and AW-1 were treated as minimum age estimates, whereas PD40 and the combined U-series-ESR age on WU1601 were treated as finite age estimates. Since the Al Wusta sequence accumulated over a short period of time, and contains only five finite ages (and three minimum ages), the General Outlier Model31 was unable to function, and instead a simpler model using agreement indices was employed. This analysis yielded Amodel (76) and Aoverall (79) values well in excess of the generally accepted threshold (6031), with only one age yielding an individual agreement index below this threshold (PD17, 51). These data indicate that no ages should be excluded from the model, and that the age model itself is robust. The Bayesian sequential model yielded an age for the end of Phase 1 of 93.1 ± 2.6 ka (1 σ uncertainties), while Phase 2 yielded start and end dates of 92.2 ± 2.6 ka and 90.4 ± 3.9 ka respectively. The end date for phase 2 should be treated as a maximum value since no overlying material is present, precluding the possibility of further constraining the end of this phase.

Stratigraphy and sedimentology

Sediment analysis

Bulk samples (in the form of coherent blocks) were taken at 10 cm intervals through each of the marl beds in four sections (Fig. 1C and Supplementary Figures 13 and 14). Each block was air-dried and subsamples (ca 0.5 g) were removed, powdered and analysed for percentage carbonate content using Bascomb calcimetry, which measures the volume of carbon dioxide liberated from a known sample mass during reaction with 10% HCl65. Thin sections were prepared from fresh sediment blocks. The sediments did not require acetone treatment as they were already dry and, due to their permeability, were impregnated with a bonding resin. Standard thin section preparation was then carried out using techniques developed in the Centre for Micromorphology at Royal Holloway, University of London66. Thin sections were analysed using an Olympus BX-50 microscope with magnifications from 20x to 200x and photomicrographs were captured with a Pixera Penguin 600es camera. A point-count approach was used to produce semi-quantified data from the thin sections, based on counting micro-features at 3 mm intervals along linear transects 1 cm apart. Kemp67, Stoops68 and Alonso-Zarza42 were referred to when identifying features. X-ray diffraction analysis (XRD) was carried out in the Department of Earth Sciences (Royal Holloway, University of London). Powdered samples were analysed on a Philips PW1830/3020 spectrometer with copper Kα X-rays. Mineral peaks were identified manually from the ICDD Powder Diffraction File (PDF) database. The methods and results are described further in Supplementary Information 5.

Diatoms

Sample preparation. Samples were analysed using the standard method of Renberg69 (Supplementary Information 5). Thus, all samples were treated with 30% H2O2 and 5% HCl to digest organic material and remove calcium carbonate. Distilled water was added to dilute the samples after heating, which were then stored in the refrigerator for four days to minimise further chemical reactions. The samples were rinsed daily and allowed to settle overnight. A known volume of microspheres was added to the supernatant after the last rinse to enable calculation of the diatom concentration70. The slides were air-dried at room temperature in a dust free environment before mounting with Naphrax diatom mountant. Diatom taxonomy followed Krammer and Lange-Bertalot71–73 and taxonomic revisions 74,75 with at least 300 valves enumerated for a representative sample at x1000 magnification.

Numerical analysis. Prevalent trends in the diatom assemblage were explored using ordination analyses using CANOCO 4.5 of ter Braak and Šmilauer76. Detrended Correspondence Analysis (DCA77) with detrending by segments and down-weighting of rare species was used to investigate taxonomic variations within each site and to determine whether linear or unimodal models should be used for further analyses. If the gradient length of the first axis is <1.5 SD units, linear methods (Principle Component Analysis, PCA) should be used; however, if the gradient length is >1.5 SD units, unimodal methods (Correspondence Analysis) should be used78. Detrended Canonical Correspondence Analysis (DCCA79) was also used to show changes in compositional turnover scaled in SD units. Therefore, variations in the down-core DCCA first axis sample scores show an estimate of the compositional change between samples along an environmental or temporal gradient. Depth was used as the sole constraint as the samples in each site are in a known temporal order80. The dataset was square-root transformed to normalise the distribution prior to analyses. Optimal sum-of-squares partitioning81 with the program ZONE82 and comparison of the zones with the Broken-stick model using the program BSTICK83 were used to determine significant zones. The planktonic: benthic ratio, habitat summary, concentration and the F index (a dissolution index84) were calculated for all the samples.

Stable isotopes

It is common practice, when analysing the δ18O and δ13C values of lacustrine/palustrine carbonates to either: 1) sieve the sediment and analyse the <63µm fraction, or 2) use the microstructure of the sample, as identified under thin section, to identify pure, unaltered fabrics, which can then be drilled out and analysed85. The former procedure ensures that the analysed fraction comprises pure authigenic marl (rather than a mixture of osctracod, mollusc, chara and marl components that will contain different isotopic values). The latter is done to ensure that any carbonate that has been affected by diagenesis is sampled. Neither of these approaches were carried out here as; 1) microfabric analysis showed no evidence for diagenesis (although some of the samples are cemented the cement makes a negligible component of sample mass), and 2) some of the samples have incipient cementation, which means that they cannot be sieved. Bulk carbonate powders were consequently analysed for δ18O and δ13C. To show that the analysis of bulk samples had no impact on the derived isotopic data, samples that were friable enough to be sieved were treated with sodium hexametaphosphate to disaggregate them and then homogenised and separated into two subsamples for isotopic analysis; (1) a sieved <63µm fraction and (2) a homogenised bulk sample. The resulting isotopic data showed no difference between the δ18O and δ13C values of the sieved and bulk samples (Supplementary Figure 13b), highlighting that the homogenous and unaltered nature of the material results in bulk carbonate isotopic analysis generating valid data. Two samples were taken from different locations of each sampled block to generate a larger dataset of independent samples. The δ18O and δ13C values of each samples were determined by analysing CO2 liberated from the reaction of the sample with phosphoric acid at 90°C using a VG PRISM series 2 mass spectrometer in the Earth Sciences Department at Royal Holloway. Internal (RHBNC) and external (NBS19, LSVEC) standards were run every 4 and 18 samples respectively. 1σ uncertainties are 0.04‰ (δ18O) and 0.02‰ (δ13C). All isotope data presented in this study are quoted against the Vienna Pee Dee Belemnite (VPDB) standard.

Vertebrate fossil analyses

Each fossil specimen was identified to lowest taxonomic and anatomical level possible (Supplementary Figure 20, Supplementary Table 19 and Supplementary Information 6). Taxonomic identification and skeletal element portions were determined based on anatomical landmarks, and facilitated by comparisons with the Australian National University Archaeology and Natural History reference collection (Canberra), unregistered biological collections held at the University of New South Wales (Sydney), and the large mammal collections of the Zoologische Staatssammlung München (Munich). Each specimen was assigned a size category (small, medium, and large) following Dominguez-Rodrigo and colleagues86, and corresponding to the five size classes described in Bunn87, where small, medium and large denote size classes 1-2, 3A-3B and 4-6, respectively. Element abundance is reported as Number of Identified Specimens (NISP).

Each specimen was examined for modification by eye and hand-lens (10x) under both natural and high-incidence light, and examined at different angles to assist identification of fine-scale surface modifications. Where required, further examination and photography was carried out using a digital microscope (Model: Dino-lite, AM7013MZ). Morphometric data (length, breadth and width) was measured using digital callipers (Model: Mitutoyo Corp, CD-8”PMX), and specimen weights using a digital scale. Bone surface modifications were identified and recorded following standard methodologies: butchery and tooth marks88–94, burning95–96, rodent gnawing97,98, weathering99 and trampling100. Carnivore damage was categorized as pit, score, furrow or puncture, and the location noted94. Tooth mark morphometric data – short and long axes – was also recorded. Any additional modifications, i.e. polish, manganese staining, and root etching, were also reported and described. Bone breakage was recorded as green, dry, or both, following Villa and Mahieu101. Long bone circumference completeness was recorded using the three categories described by Bunn102: type 1 (<1/2), type 2 (>1/2 but < complete) and type 3 (complete).

Lithic analysis

Lithics were systematically collected during pedestrian transects and excavations of Al Wusta. This produced a total studied assemblage of 380 lithics (Supplementary Information 7). Further lithics extended for a considerable distance to the north, seeming to track the outlines of the palaeolake, but we only conducted detailed analysis on lithics from the southern part of the site, close to AW-1 and the sedimentary ridge on which it was found (i.e. south of the Holocene playa). These were analysed using the methodology described in Scerri and colleages25,103,104 and Groucutt and colleagues45,105. As well as qualitative analysis of technological features indicating particular techniques and methods of reduction, a variety of quantitative features such as dimensions, the number of scars and % of cortex were recorded. Informative examples were selected for photography and illustration. This approach allows both a characterisation and description of the assemblage and broad comparison with other assemblages from surrounding regions.

Supplementary Material

Supplementary Information is available in the online version of the paper.