Longstanding behavioural stability in West Africa extends to the Middle Pleistocene at Bargny, coastal Senegal

Middle Stone Age (MSA) technologies first appear in the archaeological records of northern, eastern and southern Africa during the Middle Pleistocene epoch. The absence of MSA sites from West Africa limits evaluation of shared behaviours across the continent during the late Middle Pleistocene and the diversity of subsequent regionalized trajectories. Here we present evidence for the late Middle Pleistocene MSA occupation of the West African littoral at Bargny, Senegal, dating to 150 thousand years ago. Palaeoecological evidence suggests that Bargny was a hydrological refugium during the MSA occupation, supporting estuarine conditions during Middle Pleistocene arid phases. The stone tool technology at Bargny presents characteristics widely shared across Africa in the late Middle Pleistocene but which remain uniquely stable in West Africa to the onset of the Holocene. We explore how the persistent habitability of West African environments, including mangroves, contributes to distinctly West African trajectories of behavioural stability.


Introduction
A total of five samples from the Bargny 1, Senegal underwent luminescence dating the Sheffield Luminescence Laboratory.
Naturally occurring potassium (K), thorium (Th), uranium (U) are the main contributors of dose to sedimentary quartz. The concentrations of these elements in the sampled sandy sediments (excluding limestone clasts) were determined by inductively coupled plasma mass spectrometry (ICP) at SGS laboratories Ontario Canada (SI Table 1.1). Elemental concentrations were converted to annual dose rates using data from Guerin and colleagues 1 . Calculations took into account attenuation factors relating to sediment grain sizes used, density and palaeomoisture (SI Table 1.1). Attenuation of dose by moisture used present-day values with a ± 3 % error to incorporate fluctuations through time (SI Table 1.1). Contributions of gamma-dose from adjacent units was modelled using data from Aitken 2 . The contribution to dose rates from cosmic sources was calculated using the expression published in Prescott and Hutton 3 (SI Table 1.1). ** K below ICP detection limits of 0.1% so 0.05% assumed. *** Dose rates adjusted for gamma dose received from adjacent units using data from Aitken (1989).

SI
As calculated the calculated dose rates are based on analyses of the sediment sampled at the present day. This assumption is only valid if no movement and/or re-precipitation of the four key elements has taken place since sediment burial. The impacts of possible post-depositional changes to carbonate content (See Nathan and Mauz 6 ) was considered. We undertook dose rate analysis both based on sediments as they were sampled and on the sediments once thy had been stripped of carbonate using concentrate HCL. This allowed three models to be developed. Firstly a model where carbonate deposition was at a similar time to sediment deposition. Secondly a linear uptake model of carbonate through burial time and thirdly a model with recent (modern) carbonate deposition. For sample Shfd19082 this gave ages of (1) 93±4 ka, (2) 108±5 ka and (3) 117±6 ka for the respective models. Despite changing the dose rate by up to 30%, no model resolved the age inversion seen within Unit 4 between samples Shfd19081 (134±5 ka) and Shfd19082 (93±4 ka). Whilst it might be tempting to report the results of model 3 for Shfd19082, thereby minimising the age reversal, having looked in detail at the sediments, whilst the % carbonate is relatively high for this sample, it is in clast form indicative of erosion and deposition of local limestone rather than post-depositional pedogenic or groundwater calcrete carbonates. Thus we had no reason to go with either model 2 or 3 and ages are reported based on the dose rate of model 1. The presence of large limestone clasts within Unit 4 could have impacted on dose-rate in another way. As dose-rates were based on ICP analysis of bulk sediment an infinity field was assumed. Dose hot/cold spots from large clasts not included in the bulk sample and ICP analysis could have affected gamma dose rates for Unit 4. Whilst in situ field measurements of gamma dose would have incorporated this aspect into any dose-rate calculations, such measurements were logistically impossible at the time of sampling. Whilst unquantified, this is unlikely to have caused a the >40% change in dose-rate required to resolve the age reversal between samples Shfd19081 and Shfd19082. Future in-situ gamma-dose measurements might resolve whether the age reversal is in part due to dose-rate over-estimation from the unsampled clastic material. Future higher resolution dating sampling and Bayesian modelling of ages might also help resolve the age reversal and help understand its cause.
The luminescence samples were prepared under subdued red lighting following the procedure to extract and clean quartz outlined in Bateman and Catt 4 . Material for dating was taken from prepared quartz isolated to a size range of 90-180 µm as the original sediment was fined grained. As this size range is sub-optimal for single grain measurement, initial measurements were made at the single aliquot level. For this, grains were mounted as a 8 mm diameter monolayer on 9.6 mm diameter stainless still disks using silkospray. The samples underwent measurement using a Risø DA-20 luminescence reader with radiation doses administered using a calibrated 90strontium beta source. Stimulation was with blue/green LEDs and luminescence detection was through a Hoya U-340 filter. Samples were analysed using the single aliquot regenerative (SAR) approach 5,6 , in which an interpolative growth curve is constructed using data derived from repeated measurements of a single aliquot which has been given various laboratory irradiations (Figure 1.1b). Five regeneration points were used to characterise growth curves, with the first regeneration point being identical to the last in order to check if sensitivity changes caused by repeated measurement of the same grains are correctly monitored and corrected for by the SAR protocol (known as the "recycling ratio"). The most appropriate preheat temperature for the samples was selected using a dose recovery preheat plateau test ( Figure S1.2). This resulted in selection of preheat temperatures of 260 °C for 10 seconds which was applied to prior to each OSL measurement to remove unstable signal generated by laboratory irradiation. De values from individual aliquots were only accepted if they exhibited an OSL signal measurable above background, good growth with dose, recycling values within ± 10 % of unity, and the error on the test dose used within the SAR protocol was less than 20 %. The samples possessed good luminescence characteristics with a rapid decay of OSL with stimulation and OSL signals dominated by a fast component (e.g. Figure S1.1a). Within the SAR protocol results which grew well with laboratory dose (e.g. Figure 1.1b). Results of different preheat temperatures (Data presented as mean with error bars based on n=3 per temperature) in recovering a ~20 Gy beta radiation dose from sample Shfd19081 (a) Given to recovered dose ratio at different preheat temperatures. (b) recycling ratio (ratio between the first and last dose point) at the different preheat temperatures. Data points in both plots are the averages of three measurements performed for each preheat temperature. Yellow shading indicates acceptable limits for dose recovery and recycling.

SI
The effects of incomplete bleaching of the sediment during the last period of transport or exposure in situ can be profound. Typically, poorly bleached sediments retain a significant level of residual signal from previous phases of sedimentary cycling, leading to inherent inaccuracies in the calculation of a palaeodose value. By plotting the replicate De data for the samples (Figure 3) some assessment of whether older or younger material has been included in the sample measurements can be made. In principle a well-bleached sample that has not been subjected to post-depositional disturbance should have replicate De data which is normally distributed and highly reproducible (see 7 : Figure S1.3; 8 ).
Where post-depositional disturbance or incomplete bleaching prior to sample burial has occurred skewing of this distribution may occur and/or replicate reproducibility may be lower 8,9 . In the case of poorly bleached material skewing should be evident with a high De tail 10 . High De tails may also be indicative of saturated samples and interpolation of the De values from the upper, low gradient part of the growth curve 11 .
As Figure S1. 3  Ages are quoted in years from the present day (2019) and are presented with one sigma confidence intervals which incorporate systematic uncertainties with the dosimetry data, uncertainties with the palaeomoisture content and errors associated with the De determination. Data shows that the samples had high reproducibility in terms of their palaeodose and therefore ages should reflect true burial ages.

Summary
This supplement describes the laboratory preparation, analysis, and quantitative analysis of plant microfossils recovered from the Bargny 1 deposits. In the field, samples were collected at 5 cm intervals across the entire profile. Because the deposits do not exhibit a uniform age-depth profile and analysis of phytoliths/pollen is labor intensive, subsampling efforts focused on evaluating changes across the MSA archaeological strata in detail (Units 4-6) and general characterization of the major units making up the overlying deposits (Units 1-3). The 31 samples processed and analyzed for plant microfossils are listed in SI Table 2.1.

Laboratory methods
Sediment samples were ground, passed through a 250-micron sieve, and placed in a shaker overnight with Calgon solution (sodium hexametaphosphate) before having sands and clays separated by settling and centrifugation-decant cycles at 277 g for 2 minutes. At this point, samples were spiked with Lycopodium spores and treated with 10% HCl in a 40º bath for 10 minutes. After centrifugationdecant cycles until pH neutral, the samples were separated by density using a solution of zinc bromide and 5% HCl with a specific gravity of 2.3 g/ml. The resulting residue was extracted in ethanol and transferred to glycerol for analysis.
This set of laboratory methods for extracting phytoliths is based on published methods for digesting terrestrial sediments and recovering pollen and phytolith microfossils 13,14 , but omits the use of a strong oxidizing agent such as peroxide or nitric acid. Judging from the well-oxidized state of the sediments, further degradation of the organic remains would only remove valuable information represented by organic microfossils (pollen microfossils, amoebas, etc.), but also permits the use of Lycopodium spores to track laboratory errors and to calculate microfossil concentrations.

Phytolith Types and Images
Phytolith and pollen microfossils were identified using a binocular light microscope at 400x-1000x magnification. Phytolith nomenclature and categories follow the International Code for Phytolith Nomenclature 15 , but we tried specifically to create sample categories consistent with Bremond and colleagues 16,17 assessment of phytoliths from surface samples across West Africa (Table 2). Samples were analyzed until at least 200 diagnostic phytoliths (Figures) were encountered (SI Table 2.2) or until 300 Lycopodium tracers were encountered. Non-diagnostic phytolith counts were tracked, but are not considered in the phytolith sum.

Plant Microfossil Recovery, Preservation, and Concentration
Concentrations of phytoliths and pollen were established by tracking Lycopodium spores encountered during analysis. Samples where 200 Lycopodium were encountered before 50 identifiable phytoliths were considered unreliable. Although the number of diagnostic pollen sampled does reach 100 in some samples, the number of damaged, fragmented, and indeterminate types is too high to pursue reliable statistical analysis of the pollen data. In general, the patterns in microfossil deposition/preservation (SI Figure 2.10) closely follow the site's sedimentology, with the highest concentrations occurring in Units 3, 4, and 5 in tandem with higher values for organic matter and sediment sorting. The concentrations and preservation of plant microfossils signals a major change in the depositional environment between Unit 3 and Unit 4.

Phytolith Results
Phytoliths from the Bargny 1 samples show a strong representation of Chloridoideae and Panicoideae phytolith types (BIL and CRO, respectively) as well as types belonging to undifferentiated monocots (ACU_BUL, BUL_FLA, and RON) (SI Figure 2.11). The frequency of SPH_ECH phytoliths is high in some of the lowermost samples, but we must use caution in interpreting this morphotype as it overlaps with sponge spicules common in coastal zones 15 . While this type was tracked during analysis, this type is excluded from the statistical analyses. Samples with excessive SPH_PSI and SPH_ECH tend to be the most weathered and have the poorest preservation. These are excluded from further analysis.

SI Figure 2.11:
Phytolith results as number of identifications per sample.

Pollen Results
Identifiable pollen was preserved in the samples and provides some insights into the kinds of vegetation cover that grew at the site during the major depositional phases (SI Figure 2.4). The lowermost samples are characterized by smaller pollen counts and low concentrations which steadily rise up to about 250 cm, where the samples yielded high pollen counts and a greater concentration of pollen. For much of the middle section of the deposits, the pollen are often damaged and although identified types are high, these tend to be types that are easily distinguished (Poaceae/Amaranthaceae). There's a modest improvement in pollen preservation/concentration between 60-90 cm as well. What we do see in these samples are examples of types introduced by trade-wind activity (Olea-type & Pinaceae), riparian/estuarine trees/shrubs (Avicennia, Syzygium), and regional Sudanian woodland pollen types (Acacia, Celtis, Trema). By identifying components of the local vegetation and broader atmospheric processes, the palynological results provide some important context for the phytolith results.

Summary Plots by Taxonomic Affiliation/Vegetation Cover
Assigning the phytolith results to taxonomic groups provides some limited insights (SI Figure 2.13). Grass phytoliths are dominant in the Bargny1 phytolith record and follow the overall trends seen in grass pollen. Panicoid grasses are somewhat more common in the lowermost sedimentary units and become rarer across the upper three units, while Chloridoid grasses are generally more common with peaks in Unit 4 and one sample in Unit 2. Although the sphereoid-echinate category might be artificially inflated by similar-looking sponge fossils, it could also represent more palm cover in Units 5 and 6. The representation of different pollen types primarily reflects preservation, but it shows two peaks in grass influx in Unit 5 and Unit 3 followed by peaks in xeric pollen types (predominantly Amaranthaceae) in Unit 4 and Unit 2. Trade wind activity is evident in Unit 6 and Unit 3, while wetland herbs and riparian shrubs/trees are most common in Units 5 and 6. Wetland herbs reappear at low frequencies in Unit 4, Unit 2, and Unit 1. Trees from the broader Sudaneo-Guinean zone are present in most of the samples. Figure 2.13: Phytolith and pollen percent by taxonomic affiliation and vegetation type, respectively.

Quantitative Comparisons with Surface Samples
Studies of surface soil phytoliths conducted by Bremond and colleagues 16,17 provide an important resource for evaluating archaeological phytolith assemblages. We use Principal Component Analysis (SI Figure 2.14) and Minimum Square-Chord Distances (SI Figure 2.15) to compare the phytolith spectra and to gain insights into the range of environments represented at Bargny 1. PCA was chosen because of the "predict" functionality in R, which is not available for other types of factor analysis (constrained correspondence, etc.). To look for changes in rainfall or the boundaries of major vegetation zones, we applied MSCD to establish the most similar sets of modern samples and plot their relevant climatic (rainfall) or geographic (latitude) representation across the depositional sequence. Comparing the Bargny 1 phytolith spectra with Bremond and colleagues 16,17 using MSCD produced some limited insights. Although the best-match approach shows some major downshifts in moisture (SI Figure 2.15) and reconstructed latitude in Unit 5 and Unit 4 (black dots with solid lines), the range of similar samples (yellow to red circles) is especially broad in the lowermost samples. In tandem with the presence of pollen from estuarine environments, which are not represented in Bremond and colleagues 16,17 surface sample studies, this signal may be interpreted two ways. First, riparian conditions may be bringing phytoliths from a larger source area to the site. Second, phytoliths primarily distinguish between two grass subfamilies in terrestrial settings, but the correlations between environmental conditions and the dominance of Chloridoid or Panicoid grasses does not hold up in riparian settings.

Phytolith Indices for Savanna Types and Grass Stress
We can also compare the values for indices used by Bremond and colleagues 16,17 to evaluate aridity/savanna type (Iph) and grass water stress (Fs) (SI Figure 2.8). These values do a good job of discriminating the boundaries between the Sudanian and Guinean zones (Iph values) and the Sahelian and Saharan zones (Fs). These values allow us to look for potential turnover between major vegetation formations.

Synthetic Results
To provide an overview of the scale and direction of past vegetation change represented in the Bargny 1 deposits, we plotted the results for four different indices used above, grouping their results by stratum (SI Figure 2.17.). The indices include Fs, Iph, Grass-Amaranthaceae pollen, and wetland/riparian/mangrove pollen. As expected, the Grass/Amaranthaceae ratio is not especially informative, but the other indices show a clear turnover after Unit 4, when the depositional environment shifts from near-coastal estuarine conditions to halophytic dry coastal plain. We further synthesize the results (SI Figure 2.18) by plotting indices for water stress, short vs. long grass savanna, and riparian/estuarine pollen influx by stratum showing both the distribution of sample results (points) and their collective behavior (kernel density plots). In Unit 5 and Unit 6, riparian/estuarine conditions are associated with a distinct compression of variability in Iph and Fst values towards tall grass and low water stress conditions. By Unit 3 there is a turnover in the depositional environment and the phytolith signal shows more influence from short-grass savanna and higher water stress. This does not necessarily imply that climatic conditions are drier in Unit 3 compared with Unit 6. Given the chronology of the deposits, the riparian-influenced zone falls sometime during MIS-6, when regional conditions are arid 18,19 and jet-stream activity is high. The presence of Pinaceae pollen (an indicator of jet stream activity) in Bargny's MIS-6 deposits is consistent with this regional picture. Thus, it is possible that local hydrogeological factors are obscuring other signals of aridity during MIS-6 at Bargny. Sometime before the Last Glacial Maximum there is a turnover in the major depositional system at Bargny, although some local spring activity is still visible in the persistent presence of Typha.

SI Figure 2.18:
Synthetic results showing wetland, water stress, and precipitation proxies with kernel density plots for each stratum.

Regional Comparison and Context
Bargny currently lies at the margin of the Sahelian and Sudanian vegetation zones. Sahelian vegetation extends from 15º N -18º N along the Atlantic coastline and is characterized by wooded grasslands and deciduous bushland, dominated by short-grasses from the Chloridoideae subfamily. Phytolith composition from surface samples shows greater contributions of flabellate and saddle-shaped morphotypes, reflecting the greater density of short-grasses and higher levels of water stress 16,17 . Sudanian woodlands occupy a narrow band of latitudes along the coastline near 15º N, but the latitudinal distribution progressively increases to below 10º N further inland. Annual grasses from the Panicoideae are dominant in the tall grass savanna communities in this region. Flabellate and saddleshaped phytoliths are less common in surface samples from these environments 16 .
The woody component of the Sahelian zone includes Acacia and Mitracarpus, both of which are common in the pollen spectra of surface samples from the region 19 . Woody vegetation makes a greater contribution to Sudanian vegetation cover. Typical woody trees include Adansonia and trees in the Combretaceae (e.g. Terminalia), but many species are shared with the Guinean rain forest zone which extends south along the Atlantic coast from 14º N. Common palynological signals of Sudanian-Guinean woody types are Celtis, Elaeis guineensis, Alchornea, and Dodonaea 19,20 . The grassy components of the Sudanian woodland and Guinean forest zones is also similar, generally showing low frequencies of flabellate and saddle-shaped phytoliths 16 .
Studies of pollen 19 and phytoliths 16 There are no accepted direct dates linking the Bargny 1 deposits with MIS-5 or MIS-4. Unit 4 caps the archaeological deposits and its phytolith and pollen signals both point to drier conditions. Although the water stress index (Fst) is low, there is turnover to drier grasses dominated by Chloridoideae and identifiable pollen shows a stronger terrestrial signal dominated by Amaranthaceae-Chenopodiaceae. Low North Atlantic SSTs and intense trade wind circulation 18 contribute to a pronounced dry trend through MIS-4 (74-60 kya). This phase is missing from some of the benchmark marine pollen records from West Africa 23 , but other records show maximum extension of C4 carbon sources 22 and low influx of pollen representing humid vegetation types 24 .
The Last Glacial Maximum (MIS-2) is characterized by cooler NA SSTs and reduced extension of the ITCZ 27 from 24-11.6 kya. Dupont and Agwu 23 find a major extension of "grass-rich dry open forest" (ibid:169) at this time and Hooghiemstra and colleagues 19 find evidence for expansion of the Sahara at this time. Unit 3 broadly follows the trajectory of MIS-3 to MIS-2, moving from intermediate to maximum aridity by the upper samples, which may overlap with the Last Glacial Maximum. The turnover to warm NA SSTs and stronger monsoonal precipitation after the Younger Dryas (ca. 11.6 kya) is widely recognized as a major humid event in terrestrial records and marine records alike 20,27 . However, Dupont and Agwu 23 and Castaeñeda and colleagues 22 both find a more limited expansion of forest pollen and sources of C3 carbon in marine records at this time compared with prior warm/humid events in MIS-3 or MIS-5. Holocene conditions at Bargny are captured by two samples, showing a return to less arid conditions during the early Holocene (Unit 2) and rather dry conditions in the near-surface sediments.

Supplementary Information 3: Archaeological Comparisons and Refugia Analysis
Data used in comparisons between Bargny and other MSA assemblages across Africa dated to MIS 6 are shown in Table SI3.1, with data for comparisons between Barny and other dated MSA assemblages from West Africa shown in Table SI3.2 Blinkhorn and colleagues 28 have recently presented a model for the distribution of refugia for human populations in Africa during the Late Pleistocene. This has highlighted the presence of a region of West Africa in which modelled Late Pleistocene precipitation 29 has remained consistently within the 68% confidence interval of mean annual precipitation (248-1403 mm) associated with mobile huntergatherer population distributions 30 . Here, we repeat this analysis, focusing on the timeframe in which Middle Stone Age assemblages are presently known to occur in West Africa (160-11ka), and concentrating on the 68% confidence interval proposed by Blinkhorn and colleagues 28 , which effectively encapsulated Late Pleistocene occupations of eastern Africa, as well as more broadly reflecting longer term MSA distribution patterns in the region during the Middle and Late Pleistocene, as identified by Timbrell and colleagues 31 .
Overall, this analysis spanning 160-11ka chronologically extends the trends observed for the Late Pleistocene, highlighting a discrete region in West Africa that is persistently within the habitable precipitation bracket throughout this extended timeframe, termed here the Senegambian refugia (SI Figure 3.1). Whilst the overall size of the Senegambian refugia is comparable to that previously identified for the Late Pleistocene, E-W connectivity with eastern Africa is notably more limited. A N-S cline in precipitation is observed in the Senegambian refugia, with increased precipitation observed in the south, whereas E-W clines are observed for both standard deviation and coefficient of variation of precipitation is seen, suggesting greatest stability in precipitation occurs in the western part of the refugia. Tropical Xerophytic woodland is the most commonly occurring habitat type, with greater stability of open habitat types present in the northern half of the Senegambian refugia contrasting with higher amplitude of habitat change and diversity seen in the southern half.     (160-11ka). These support the persistent habitability of the Senegambian refugia from the initial occupation at Bargny 1 to the onset of the Holocene and youngest MSA occupations at Saxomununya, and offer insight into the role of environmental stability on observed patterns of technological homogeneity.