First evidence of an extensive Acheulean large cutting tool accumulation in Europe from Porto Maior (Galicia, Spain)

We describe a European Acheulean site characterised by an extensive accumulation of large cutting tools (LCT). This type of Lower Paleolithic assemblage, with dense LCT accumulations, has only been found on the African continent and in the Near East until now. The identification of a site with large accumulations of LCTs favours the hypothesis of an African origin for the Acheulean of Southwest Europe. The lithic tool-bearing deposits date back to 293–205 thousand years ago. Our chronological findings confirm temporal overlap between sites with clear “African” Acheulean affinities and Early Middle Paleolithic sites found elsewhere in the region. These complex technological patterns could be consistent with the potential coexistence of different human species in south-western Europe during the Middle Pleistocene.


SI Geomorphologic and stratigraphic context
The Miño River is the most important hydrographic network in NW Iberia with an area of ~17,000 km 2 and a length of 350 km (Fig. S1). Within the lower Miño River basin, and in the surroundings of Porto Maior site, 9 stepped fluvial terraces have been identified, with the following (relative) heights above river level estimated during the summer: T1 (+4-7 m), T2 (+13-17 m), T3 (+21-29 m), T4 (+30-39 m), T5 (+45-51 m), T6 (+53-61 m), T7 (+65-77 m), T8 (+78-89 m) and T9 (+91-108 m) (Fig. 1, Fig. S2). These terraces were formed due to the effects of tectonic uplift and climatic fluctuations during the Plio-Pleistocene. The sedimentary deposits of these terraces are characterized by thick sequences made of clast-supported crudely-bedded gravels, mainly quartzite and quartzes with minor components of granites (Gh facies), with internal planar crossbedded facies (Gp) or cross-bedding (Gt), that alternate with fine to very coarse sands or pebbles with planar (Sp) or cross-bedded (Sp) structures. Also recognized are finecoarse Sm facies with massive structure or faint laminations, similar to Fsm facies of silts and muds with massive structure, occasionally displaying mud cracks (Fm facies). The facies architecture indicates that the fluvial style is of a braided river system, with deposits dominated by gravels 1,2 .  In total, 5 stratigraphic levels have been identified at Porto Maior, named PM1 to PM5 from bottom to top. PM1 and PM2 are clast-supported gravels and most of the sand matrix fraction is coarse to very coarse (0.5-2mm) (Fig. S3). The clay fraction amounts to between 20 and 25 %, and has originated from processes of illuviation (horizon Bt) that developed after the sedimentation of level PM3. The colours of the sand-clay matrix are yellow (10YR 5/6) or reddish brownish (5YR 5/4). Level PM2 is more affected by reduction processes, and exhibits a pale yellow (2.5Y 7/4) and olive (5Y 5/6) colour. PM2 is a gravel facies Gh or Gp with horizontal planar stratification and clast imbrications.
Level PM3 has a uniform textural composition, and is composed of muds with percentages of sand that range between 35 % and 40 %; although fine to median sand (between 0.5 and 0.05 mm) dominates. In these massive muddy deposits (Fsm and Fm facies), levels of low-thickness gravel are recognized associated with small flat channels, of limited lateral development. These channels would have traversed the overbank facies deposits, which in turn would have been affected by waterlogging and argilluviation  were collected from a perpendicular outcrop, a few meters away from the main excavation area. Drawing and photo by M. Duval.

Methods
Sediment samples were prepared in the laboratory under conditions of limited illumination following the standard procedure at CENIEH 4 . The 100-200 µm size fraction was collected after wet sieving. HCl (36%) was used to dissolve carbonates and H2O2 (30%) to eliminate organic matter. Heavy minerals and feldspars were removed with Sodium Polytungstate solutions at d=2.72 and d=2.62 g/ml, respectively. Then, magnetic minerals were eliminated using neodymiun magnets. The resulting samples were treated with HF (40%) for 40 minutes to eliminate the remaining feldspars and to etch quartz grains. Finally, HCl (18%) was added in order to remove any soluble fluoride. ESR measurements were carried out at CENIEH (Burgos, Spain), with an EMXmicro 6/1Bruker X-band ESR spectrometer coupled to a standard rectangular ER 4102ST cavity. To ensure constant experimental conditions over time, the temperature of the water circulating in the magnet is controlled and stabilized at 18 ºC by a water-cooled Thermo Scientific NESLAB ThermoFlex 3500 chiller, and the temperature of the room is kept constant at 20 ºC by an air conditioning unit. ESR measurements were performed at low temperature (~90 K) using a ER4141VT Digital Temperature control system based on liquid nitrogen cooling. Further details about the setup and about its stability over time can be found 5 .
In accordance with the Multiple Centre method defined by 6 , the ESR signals of both the Al and Ti centres were measured. For the first one, the following acquisition parameters were used: 10 mW microwave power, 1024 points resolution, 20 mT sweep width, 100 kHz modulation frequency, 0.1 mT modulation amplitude, 40 ms conversion time, 10 ms time constant and 1 scan. In contrast, the ESR signal associated to Ti centre was measured as follows: 5 mW microwave power, 1024 points resolution, 20 mT sweep width, 100 kHz modulation frequency, 0.1 mT modulation amplitude, 60 ms conversion time,10 ms time constant and 2 to 30 scans (depending on the aliquot). Each of the 14 aliquots (one natural, one optically bleached and eleven gamma irradiated aliquots) of a given sample were measured 3 times after a ~120° rotation in the cavity for both Al and Ti signals in order to consider angular dependence of the signal due to sample heterogeneity. The only exception was MIN1401: the very low ESR intensity of the Ti signals required a higher number of scans (up to 30), resulting in an already very long measurement time (>5 hrs). Consequently, no rotation were considered for this sample.
Then, all measurements were repeated three times over distinct days in order to check the repeatability of the DE values. Consequently, three sets of 13 data points were obtained for each signal measured (Ti and Al) in a given quartz sample.
The ESR intensity of the Al signal was extracted from peak-to-peak amplitude measurements between the top of the first peak (g=2.0185) and the bottom of the 16th peak (g=1.9928) 7 . Following the conclusions from 8 , the ESR intensity of the Ti-Li centre was preferentially evaluated by measuring the peak-to-baseline amplitude around g=1.913-1.915 (option D). Peak-to-peak amplitude measurement between g=1.979 and the bottom of the peak at g=1.913 was also performed for comparison (option A). The intensity of the Ti-H centre was evaluated by taking the peak-to-baseline amplitude around g=1.915.
For each aliquot, ESR intensities of Al and Ti centres were corrected by the corresponding receiver gain value, number of scans, mass and a temperature correction factor 5 . The fitting procedures were carried out with the Microcal OriginPro 9.5 software using a Levenberg-Marquardt algorithm by chi-square minimization. For the Al centre, an exponential+linear function (EXP+LIN) was fitted through the experimental points (see equation in 4 ), and data were weighted by the inverse of the squared ESR intensity (1/I 2 ). DE values were obtained by extrapolating the EXP+LIN function to the residual intensity (so-called Total bleach method, 9 ). For the Ti centre, we used the function labelled as Ti-2 in 8 , in order to describe the non-monotonic dose dependence of the ESR signal at high doses. Data were weighted by the inverse of the squared experimental error (1/s 2 ) and DE values were obtained by back extrapolation to the Y=0.
For each sample, final dose response curves (DRCs) were obtained by pooling all the repeated measurements in a single plot, as recommended by 10 .
The total dose rate value was derived from a combination of in situ and laboratory measurements. External gamma dose rate were derived from in situ measurements by using the "threshold technique" 11 . For each dated samples, the corresponding radioelement (U, Th, K) concentrations in the sediment were determined by ICP-MS analysis of about 5g of dry raw sediment. In addition, ~150 g of this same raw sediment, previously dried and powdered, were analysed by High Resolution Gamma Spectrometry (HRGS) using a Canberra Extended Range (XTra) HpGe detector in order to identify possible disequilibrium in the U-238 decay chain. Concentration values were used to derive external alpha and beta dose rate components using the dose rate conversion factors from 12 . Dose rate values were calculated assuming a mean grain size of 150 µm, and an assumed thickness removed by HF etching of 20 µm. Internal dose rate was assumed to be 50+30 uGy/a, based on the work from 13 and assuming an alpha efficiency k of 0.15+0.10 14 . Values were corrected with beta and alpha attenuation values for spherical grains 15,16 and water attenuation formulae from 17 . Current water contents were evaluated in the laboratory by drying the sediment at 50ºC in an oven during three weeks. Results vary within relatively narrow range from 8.5 to 13.2% (wet weight) among the samples. For a matter of consistency with the luminescence dating procedures, a value of 20+5% (wet weight) was considered in the ESR age calculation of all but one sample, which is equivalent to the 20-28% (dry wet) used for the OSL and pIR-IR samples. For the uppermost sample MIN1403 (unit PM1), a water content value of 13+5% was considered instead, which is also similar to the value of 15±5% (dry weight) adopted for the luminescence sample collected from this unit (PM16-4). The cosmic dose rate was calculated using formulae from 18      In accordance with the basic principles of the Multiple Centres approach, this suggests that for these samples the signal of the Al centre has been incompletely reset during sediment transportation. Consequently, the Ti-Li centre most likely provides the best estimation for the burial dose of the samples.

Dose rate evaluation
Sediment was analysed in situ as well as in the laboratory by a range of different techniques. ICP-MS and HRGS analyses overall provide highly consistent results ( Table   S5). Some of the differences observed may be due to some inherent variability within the sediment, as ICP and HRGS analyses were performed on ~10 g and ~150 g of raw Ti-Li age ratio AW/MW 1.06 1.09 1.14 1.14 1.13 1.06

SI Luminescence dating
In total, four samples were collected for luminescence dating at the site of Porto Maior ( Fig. 1). Two samples, PM16-3 and PM16-5, were collected from fluvial sediment layers located immediately below the excavated archaeological horizon, within lithostratigraphic unit PM3. The lowermost sample, PM16-3, was collected ~1.5 m below the archaeological layer and near the base of a massive to weakly bedded sandy deposit.
Sample PM16-3 was also located ~20 cm above the contact with the underlying deposit μm quartz fraction of sample PM16-4 were prepared following the same procedure, with the exception that heavy liquid densities of 2.62 g/cm 3 and 2.72 g/cm 3 were used for separating feldspars and heavy minerals, respectively. The selected quartz grain fraction was etched with 48% hydrofluoric acid for 40 minutes. As with K-feldspars, samples were re-sieved with a 63 µm sieve after the etching procedure.

Environmental dose rate estimation
Dose rate assessments were made using a combination of in situ gamma spectrometry measurements and low-level beta counting (Table S5). Field gamma spectrometry measurements were performed at each luminescence dating sample position immediately after removal of the PVC tubes. Elemental concentrations of K, U and Th were determined from the field gamma spectra using the 'energy windows' method described in 21 . Additional sediment samples were collected from the area around each luminescence dating sample position for beta dose rate assessments (beta counting), water content evaluations, and high-resolution gamma-ray spectrometry (HRGS) measurements. Low-level beta counting was performed on dry and homogenised sediment using a Risø GM-25-2 beta counter 22 (Table S6). Cosmic-ray dose rates were calculated using the approach described in 18 . Internal dose rate contributions for K-feldspar grains have been estimated using an assumed internal 40 K content of 12.5 ± 0.5% 23 and 87 Rb content of 400 ± 100 ppm 24 .
The beta, gamma and cosmic-ray dose rates have been corrected for estimated longterm water contents of each sample 25,26 . The present-day sediment water contents ranged between 10 and 19% of dry sediment weight (Table S13) but they are not considered to be entirely representative of long-term moisture conditions at the site because the sediment profiles had been exposed for 3 years in the excavation area (samples PM16-6 and PM16-4), and for more than 50 years in the abandoned railway trench (samples PM16-3 and PM16-5); and thus they had partially dried out prior to sampling. To determine more suitable long-term sediment moisture contents, we have adopted conservative estimates based on 50% present-day saturated water contents for each luminescence sample. A 1σ relative uncertainty of 20% has been assigned to the long-term moisture estimates to accommodate any potential variations in hydrologic conditions during burial (e.g., changes in effective sediment moisture between glacial and interglacial cycles). This approach yielded long-term sediment moisture contents of 20-29% for samples PM16-3, PM16-5 and PM16-6 (Table S5), which overlap with published values used in reliable, known-age luminescence dating studies of Middle Pleistocene fluvial and alluvial deposits from Spain (e.g., [27][28][29][30]. For sample PM16-4, which comes from the uppermost aeolian unit (PM5), the water content value used to calculate the final ages is based on values adopted in other luminescence dating studies of wind-blown (loess) sequences from Europe. These values range from 10% to 20%, and so we applied a value of 15 ± 5% for samples PM16-4 [31][32][33][34][35][36][37] . Table S13 compares the final pIR-IR ages obtained using our preferred estimates of the long-term water contents (which account for recent loss of moisture due to profile exposure) and the ages that would be obtained using the present-day (as measured) water content values of each sample. It can be observed that the corresponding ages for each sample are not statistically different at 2σ, and therefore we conclude that our age estimates and interpretations are not sensitive to our preferred choice of long-term water content. The adopted long-term moisture estimates, which are expressed as percentages of dry sediment weight, are also consistent with the long-term water contents for the ESR dating samples (20 ± 5%), which are expressed as percentages of wet sediment weight.

Instrumentation and equivalent dose (De) measurements and estimation
Post-infrared infrared stimulated luminescence (pIR-IR) measurements were carried out using a Risø TL/OSL-DA-20 reader equipped with a calibrated 90  (ii) the recuperation ratio, calculated as the ratio of the sensitivity-corrected 0 Gy dose point (L0/Tx) to the sensitivity-corrected natural (Ln/Tn), was <5%; (iii) the sensitivitycorrected natural signal intercepted the sensitivity-corrected dose-response curve or it did not intercept the saturated part of the dose-response curve (i.e., the Ln/Tn value was equal to, or lower than, the Imax saturation limit of the dose-response curve at 2σ); (iv) the dose-response curve did not display anomalous properties (e.g., zero or negative responses with increasing dose) and resulted in good Monte Carlo fits.
Single-grain OSL De values were measured using the SAR protocol shown in Table S7.

pIR-IR dose recovery tests and signal characteristics
To determine the most suitable pIR-IR measurement and preheat conditions for the Porto Maior samples, we undertook dose recovery tests using the pIR-IR225 and pIR-IR290 SAR protocols shown in Table S7. Ten 160-grain K-feldspar aliquots of sample PM16-6 were  A representative pIR-IR225 decay curve and sensitivity-corrected dose-response curve is shown in Fig. S9a-b. The pIR-IR225 decay curves of these samples typically decrease by ~90% within the first 30 s of stimulation and are optimally fitted with a single saturating exponential plus linear function. All the De values were obtained from the region of the dose-response that was not in saturation when using this type of fitting function. All of the measured aliquots passed the SAR quality assurance criteria outlined in the previous section. Fig. S9c-d shows the OSL dose-response and decay curves obtained for a grain of sample PM16-4. The example shown represents a typical grain accepted after applying the SAR quality assurance criteria. The OSL signal is fast-decaying and reaches background within 0.6 s of stimulation, while the dose response curve is best fitted with an exponential function.

Residual dose assessments
In order to assess the bleaching properties of the pIR-IR225 signal for these samples, and Consequently, we do not consider the low g-values recorded here to be indicative of the need for pIR-IR age corrections. The consistency of the existing pIR-IR225 ages and the replicate single-grain OSL and quartz ESR ages at Porto Maior (Fig. 1) similarly does not support the need for additional fading corrections in this context. a Water content used for calculating environmental dose rates, expressed as % of dry mass of sample and assigned a relative uncertainty of ± 20%. For samples PM16-3, PM16-5 and PM16-6 the long-term water contents are 50% of saturated values. For sample PM16-4 the long -term water content value was set to 15 ±5%. b Radionuclide concentrations and specific activities have been converted to dose rates using the conversion factors given in 57 and 58 , making allowance for beta-dose attenuation 15,59 . c Gamma dose rates were calculated from in situ measurements made at each sample position with a NaI:Tl detector using the 'energy windows' method detailed in 21 . d Beta dose rates were calculated using a Risø GM-25-5 low-level beta counter 22 , after making allowance for beta dose attenuation due to grain-size effects and HF etching 15 . e Cosmic-ray dose rates were calculated according to 18 and assigned a relative uncertainty of ± 10%. f Assumed internal (alpha plus beta) dose rate for K-feldspar grains, the Internal alpha and beta dose rate contributions from 238 U and 232 Th were calculated using assumed concentrations of 0.15 ± 0.03 ppm and 0.35 ± 0.07 ppm, respectively, based on modal values obtained by 60 and similar values obtained by [61][62][63] . An a-value of 0.09 ± 0.03 was used to estimate the internal alpha dose rate contributions from these 238 U and 232 Th concentrations based on published estimates obtained for a wide range of k-feldspar samples -e.g. 27,64-68 -. The assumed internal (alpha plus beta) dose rate for the quartz fractions is based on published 238 U and 232 Th measurements for etched quartz grains from a range of locations 60,69-71 and an alpha efficiency factor (a-value) of 0.04 ± 0.01 64,72 . g Internal dose rate of feldspar grains arising from 40 K and 87 Rb concentrations were calculated from assumed values of 12.5 ± 0.5% 23 and 400 ± 100 ppm 24 , respectively. h Mean ± total uncertainty (68% confidence interval), calculated as the quadratic sum of the random and systematic uncertainties.  . Values shown are the specific radionuclide activities (Bq kg -1 ) and daughter-to-parent ratios.

De results and ages
Step SAR pIR-IR225 SAR pIR-IR290 Step SAR single-grain OSL Signal      55 .
Critical skewness values are taken to be equivalent to the standard error of skewness score (68% C.I.) for multi-grain aliquot De datasets and twice the standard error of skewness score (95% C.I.) for single-grain De datasets, following the results of sensitivity analyses performed by 55 and 73 . d Mean ± total uncertainty (68% confidence interval), calculated as the quadratic sum of the random and systematic uncertainties. Total uncertainty includes a systematic component of ± 2% associated with laboratory beta-source calibration.  55 . Critical skewness values are taken to be equivalent to the standard error of skewness score (68% C.I.) for multi-grain aliquot De datasets and twice the standard error of skewness score (95% C.I.) for single-grain De datasets, following the results of sensitivity analyses performed by 55 and 73 .
c Mean ± total uncertainty (68% confidence interval), calculated as the quadratic sum of the random and systematic uncertainties. Total uncertainty includes a systematic component of ± 2% associated with laboratory beta-source calibration. Table. S 13. Ages obtained for the Porto Maior luminescence dating samples using the present-day 'as measured' water content values and the assumed long-term water content values. The latter have been used to derive the final age estimates and equate to 50% of the saturated water contents for samples PM16-3, PM16-5 and PM16-6, and a mean published value of 15 ± 5% for the uppermost aeolian sample PM16-4. Measurements were made on 1500-grain aliquots containing quartz grains and using a SAR protocol (see Table S7). The test-dose preheat used for this experiment was set to 160°C for 10 s for all measurements, while the regenerative dose preheat was varied between 200°C and 260°C (for 10 s) as shown. (B) Radial plot showing the single-grain OSL dose recovery test results obtained for sample PM16-4 using a regenerative-dose preheat of 220°C for 10 s and a test dose preheat of 160°C for 10 s. The grey bar is centred on a recovery dose ratio of 1 (the weighted mean measured to given dose ratio = 0.97 ± 0.02).  Table S12 for details).

SI Site formation
Level PM4 has a thickness of ~80 cm and comprises fine sediments (22) that are indicative of a low energy sedimentation process (Fsm facies). The lithic material is located at the base of the level (Fig. 2, Fig. S12-13 Further evidence supporting the separation of the entire assemblage into sub-series comes from the size-and weight-range of the material (Fig. S15-16). The MA series is almost exclusively composed of large size pieces, while the FRA series shows a clear distribution toward smaller ranges. These differences relate to the degree of fluvial abrasion; the larger pieces are substantially less abraded than the smaller ones ( Fig.   S15-16).
Quantitative evaluation of the distribution patterns (Fig. S17), via the use of cluster analysis (k-means method 75,76 ), revealed the identification of two statistically significant groups (Fig. S17a). These groups were analysed according to the type of series. We identify group 1 as being composed of pieces of smaller dimensions (flakes, cores, waste...), and group 2 primarily made of LCTs (Fig. S17b-c). The correlation between taphonomic and spatial pattern supports the separation in two different sub-series (MA,         78 .

SI Technological features of the LCTs implements
The preferred blanks for the handaxes were pebbles (63.7 %), while the use of flakes as blanks was less frequent (25.0 %) (Fig. S16-17). The presence of indeterminate blanks is low (8.7 %). The low percentages of flake blanks is related to the favourable size and shape of the pebbles, which eliminates the need for undertaking intermediate steps when

SI Overview of extra-European sites with large LCT accumulations
Sites with large accumulations of LCTs are common in Africa and the Near East, but they were hitherto unknown in Europe ( Fig. S25; Table S17).  Fig. 3; Table S17). Commonly, these implements are associated with variable percentages of flakes, cores and waste. However, it seems that the LCTs were configured elsewhere and were subsequently introduced to the site.
Normally, the lithic assemblages are not found in direct stratigraphic relation with faunal remains 82,135,136 . Where faunal remains do exist, they appear to be coincidental accumulations rather than the product of butchery practices.
These African and Near Eastern LCT sites comprise an extensive chronological range