Hominin presence in Eurasia by at least 1.95 million years ago

Domain: earth_science
Published: 2025-01-20
Editorial outcome: major_revision

Abstract
The timing of the initial dispersal of hominins into Eurasia is unclear. Current evidence indicates hominins were present at Dmanisi, Georgia by 1.8 million years ago (Ma), but other ephemeral traces of hominins across Eurasia predate Dmanisi. However, no hominin remains have been definitively described from Europe until ~1.4 Ma. Here we present evidence of hominin activity at the site of Grăunceanu, Romania in the form of multiple cut-marked bones. Biostratigraphic and high-resolution U-Pb age estimates suggest Grăunceanu is >1.95 Ma, making this site one of the best-dated early hominin localities in Europe. Environmental reconstructions based on isotopic analyzes of horse dentition suggest Grăunceanu would have been relatively temperate and seasonal, demonstrating a wide habitat tolerance in even the earliest hominins in Eurasia. Our results, presented along with multiple other lines of evidence, point to a widespread, though perhaps intermittent, presence of hominins across Eurasia by at least 2.0 Ma.

Introduction
Current evidence for the earliest appearance of hominins outside Africa comes from the site of Dmanisi, Georgia. Dated to 1.85-1.77 Ma1, the Dmanisi assemblage includes a large number of hominin remains2, as well as lithics and evidence of hominin modification of animal remains (e.g., butchery marks3). This site clearly demonstrates a hominin presence in Southwest Asia/Eastern Europe by the Early Pleistocene (Gelasian), yet the exact timing of the initial dispersal of hominins out of Africa and the long-term success of these dispersals is debated (e.g., 4). This is especially true for Europe, where there is an ongoing discourse regarding the timing of hominin presence in the earliest appearance of hominins outside southern and northern Europe5,6. Updated chronologies for fossil localities in Asia (especially China7) and a growing number of localities that may represent fleeting or ephemeral traces of hominin activity (e.g., lithics and/or anthropogenic modifications of bones unaccompanied by hominin fossils) throughout Eurasia increasingly suggest hominins were likely present in Eurasia prior to Dmanisi (Fig. 1).

One Early Pleistocene Eurasian locality that could shed light on the initial dispersal of hominins into Eurasia is Grăunceanu, located in the Olteţ River Valley (ORV) of Romania. This region is situated in the Dacian sedimentary basin, just south of the Carpathian Mountains (Fig. 1). Deposits stem from the Tetoiu Formation, sediments of which represent multiple fluvio-lacustrine sequences that are rich in fossils8,9. This formation extends from the base of the Pleistocene to as young as ~1.3 Ma8 (see Supplementary Note 1). Grăunceanu was originally excavated in the 1960s and is one of the best known Early Pleistocene sites from East-Central Europe. Biochronological assessments indicate Grăunceanu is Late Villafranchian (~2.2–1.9 Ma) and is attributed to mammalian biostratigraphic zones MN17/MmQ19–11. At least 31 taxa (Supplementary Data 1) are identified from Grăunceanu, including mammoth, multiple species of bovids and cervids, giraffids, equids, rhinocerotids, multiple carnivore species, rodents (beaver, porcupine), ostrich, a large species of terrestrial monkey (Paradolichopithecus), and the youngest representative of pangolins in Europe11–14. Paleoecological analyses suggest Grăunceanu was a forest-steppe environment along the paleo-Olteţ river15. Other localities in the ORV include the penecontemporaneous sites of La Pietriș and Valea Roșcăi and the slightly biochronologically younger Fântâna lui Mitilan9, as well as multiple smaller localities. Though no hominin remains or in situ lithics have been identified from Grăunceanu, prior researchers have described lithics (see Supplementary Note 2 and Supplementary Fig. 1) from the nearby penecontemporaneous site of Dealul Mijlociu16, which Radulesco and Samson9 placed in a similar faunal horizon as Grăunceanu based on their biochronological comparisons.

Here we present evidence of hominin presence in Eurasia by at least 1.95 Ma in the form of cut-marked bones from the site of Grăunceanu, Romania, supported by high-precision uranium-lead (U-Pb) age estimates. We also perform high-resolution oxygen and carbon stable isotope analysis of a horse maxilla from the same site to reconstruct temperature seasonality and precipitation. We place these data into the context of the larger discussion regarding ephemeral traces of early hominin dispersals into Eurasia in the Early Pleistocene and argue in favor of a hominin presence across Eurasia by at least 2.0 Ma.

Results
Taphonomy
The total number of identified specimens (NISP) from Grăunceanu is 4983. Of these, 4524 (excluding isolated teeth and horn/antler specimens) were examined under strong, low-angled light for evidence of bone surface modifications (BSMs) by one or more observers, following protocols outlined in ref. 17. The taphonomic condition of the remains is very consistent (i.e., all bones show similar coloration and post-depositional alterations), with many bones still in articulation at the time of excavation, and the majority of specimens (>75%) preserving at least half of the original skeletal element. We interpret this assemblage as representing a single deposition, likely a low-energy overbank seasonal flood deposit. The assemblage overall shows little weathering (85.3% of specimens at weathering stage 0) or water-based alterations (only 21 specimens have polishing). These results suggest the assemblage was not subject to reworking and was not exposed for long on the surface before burial. Bone surface visibility overall is very good, with most specimens (>73%) having 75–100% of the bone surface visible. A large proportion of specimens (81.7%) present some degree of root-etching, and post-depositional damage (e.g., chipping, cracking, exfoliation, etc.) was present on 41.5% of specimens. Carnivore damage is present on 9.5% of the assemblage in the form of tooth scores, pits, crenulated break edges, or any combination of these three bone modifications.

A total of 1,189 specimens from Grăunceanu exhibit linear marks; the vast majority of these were identified as tooth (n=290), trampling (n=172), excavator (n=296), or unknown (n=411) marks (see Supplementary Notes 3-5). No evidence of percussion marks was recorded. Twenty specimens total exhibit cut marks; of these, 7 display high-confidence cut marks, 12 show probable cut marks, and 1 specimen presents both types of marks (Fig. 2, Supplementary Figs. 2–21, Supplementary Data 2); detailed descriptions of all marks can be found in Supplementary Note 5. Cut marks were identified using two methods: 1) qualitative analyses modified from 18, 19, and 2) quantitative analyses using methods outlined in ref. 20 (Supplementary Figs. 22–24, Supplementary Table 1). The eight specimens identified as having high-confidence cut marks include four tibiae, one mandible, one humerus, and two long bone fragments. All of the taxonomically identifiable specimens with high-confidence cut marks are from artiodactyls, except one small carnivoran tibia. Most specimens have two or more linear marks identified as cut marks. All high-confidence cut marks have straight trajectories, primarily transverse orientations relative to the long axis of the bone, and the same color within the mark as the main bone surface, while very few have preserved evidence of microstriations or shoulder effects. For specimens that could be identified to element, the cut marks appear in anatomical locations consistent with defleshing, especially the distal tibia. Most marks identified as high-confidence cut marks were also classified as cut marks by the quantitative analysis with high posterior probabilities (Supplementary Data 2).

The 13 specimens exhibiting probable cut marks were given lower confidence due to these specimens having more degraded surfaces (e.g., root etching, surface exfoliation), not being identified as cut marks in the quantitative analysis, and/or not presenting enough qualitative features to be definitively identified as cut marks. These marks are found on a wider variety of skeletal elements, though metapodia are the most frequent (n=5), followed by other long bones and mandibles. Probable cut marks could be a single long linear mark or a patch of up to 15 small linear marks. Most have straight trajectories, are more frequently obliquely oriented to the long axis of the bone, rarely have microstriations, and none exhibit shoulder effects.

U-Pb dating
Radiometric ages for ORV localities were estimated using laser ablation U-Pb analyses (broadly following 21 and parameters detailed in Supplementary Table 2) of dentine from seven mammal specimens from Grăunceanu, one from Valea Roșcăi, and one from Fântâna lui Mitilan (Fig. 3, Supplementary Tables 3 and 4, Supplementary Figs. 25 and 26). In Supplementary Table 5 we show both “equilibrium” age estimates, based on the U-Pb data, and also slightly younger “disequilibrium-corrected” ages, which take into consideration potential initial disequilibrium in the 238U-206Pb decay chain (see Supplementary Note 6). All subsequent discussion utilizes the latter, which provide the most robust age estimates for the Grăunceanu samples and range from 2.01±0.20 Ma to 1.87±0.16 Ma (±2σ) with an average of 1.95 Ma; ages cluster tightly around this average (Fig. 3), lending high confidence to this estimate. Specimen VRc.0001 from Valea Roșcăi returned a corrected age of 2.05±0.22 Ma, while specimen FM.0019 from Fântâna lui Mitilan returned a corrected age of 1.61±0.10 Ma. These U-Pb estimates agree with biochronological estimates for Grăunceanu and Fântâna lui Mitilan 11 and with the stratigraphic positioning of Valea Roșcăi, which has been suggested to be similar in age to Grăunceanu9 (see Supplementary Note 6). These results provide our best estimates of the minimum ages of fossil burial, as they represent final closure of the U-Pb system following post-depositional uranium exchange with the environment, which is ubiquitous in fossil teeth and bone samples. The time between deposition and closure is almost certainly variable and probably site related. These U-Pb results, coupled with the biochronology for these localities 11–14, lend high confidence to our estimated age of deposition for Grăunceanu of >1.95 Ma and potentially prior to 2.0 Ma, and indicate that Grăunceanu may be the earliest locality in Europe to show evidence for hominin activity.

Stable isotopes
We analyzed oxygen (δ18O) and carbon (δ13C) stable isotope ratios from post-weaning cheek teeth of an Equus sp. specimen from Grăunceanu (Fig. 4, Supplementary Data 3, Supplementary Note 7). We used high-resolution sampling to assess variation across the tooth crown and along the tooth row (P2-M3); this allowed us to examine seasonal fluctuations in rainfall, as equids are obligate drinkers (i.e., they depend on fresh water sources) and their diet is relatively fixed (i.e., grazing).

The δ18O values we observe in carbonate vary between 19.9‰ and 26.2‰ (Vienna Standard Mean Ocean Water; VSMOW), with an average value of 23.5‰. Curran et al. 15 reported mean δ18O values from Grăunceanu of 18.9‰ to 25.3‰ (VSMOW) and this slight difference could be explained by the fact that their values were measured on bulk samples from multiple artiodactyl specimens, incorporating enamel with an undetermined seasonal bias, or could represent variation in feeding behaviors and/or migratory habits.

We converted these δ18O isotopic values to meteoric water values using three different transfer functions 22–24. The average of annual values across the P2-P3-P4-M3 tooth series and across the three equations is −10.8±1.0‰ (VSMOW), which is lower than the present-day weighted mean at Râmnicu Vâlcea of −8.1‰. Such a large difference would result in annual temperature reconstructions of ~4°C if the equations of Rozanski et al. 25 or Pryor et al. 26 were employed, but this reconstruction is incompatible with the presence of warm-adapted fauna at Grăunceanu 13–15. Multiple factors could influence the isotopic signature of rainfall at Grăunceanu (e.g., air mass origin, rain-out history, global ice volume, moisture recycling); however, information on these factors from the Early Pleistocene is scarce and does not allow us to discuss their potential isotopic impact on regional rainfall. Thus, we argue that this lower value suggests an increased winter precipitation component to aquifer recharge compared to present. An increased winter component in the Early Pleistocene could also explain why reconstructed summer values at Grăunceanu (−7.3‰) are lower than regional values for today 27 (−3.9‰), as the horse might have consumed water with a stronger winter isotopic signal (see Supplementary Note 7).

Conversion of our enamel δ13C values to plant tissue δ13C yields values between −28.9±0.5‰ and −25.7±0.5‰ (VPDB), consistent with feeding in a woodland to woodland-mesic C3 grassland 28, and supports prior environmental reconstructions for Grăunceanu 15. The highest values are recorded in summer, which could be related to photosynthesis under drought stress 29, while the lowest values are recorded during winter (Fig. 4). Using the equation of Kohn 29 and after performing full error propagation, we obtained precipitation values between −4 and 32 mm for the driest month, and between 45–222 mm for the wettest month. This result further supports our interpretation of an increased winter precipitation component to aquifer recharge in the ORV.

Methods
Materials examined as part of this research are housed at the “Emil Racoviţă” Institute of Speleology (ISER) and the Museum of Oltenia (MO); materials from the ORV are currently split between these two institutions. These collections are curated by coauthors Petculescu (ISER) and Popescu (MO) and permissions to access the collections for this research were granted by the curators and associated institutions. Excavations at Grăunceanu were undertaken from 1960–1966, with some additional test excavations and/or collection as early as 1959 and as late as the 1970s. Nearly all excavation records and provenience data have unfortunately been lost. Publications by Bolomey 10 and Radulesco and Samson 9 describe the open-air site of Grăunceanu as a ~0.75 meter (m) thick bone bed spread over an area 90 m2, with little to no stratification of the fossil horizon. All Grăunceanu materials examined as part of this study, including those sampled for U-Pb dating and isotope analysis, come from the original excavations that took place in the 1960s and stem from this single faunal accumulation. We therefore have high confidence that these remains all represent the same depositional event. Detailed taphonomic descriptions of the ORV localities are forthcoming.

Our team has successfully reidentified the locations of multiple of the fossil localities in the ORV, and we have accessioned and/or inventoried these collections. The total number of identified specimens (NISP) for the Olteţ River Valley Assemblage is currently 5527. This includes materials from Grăunceanu (NISP=4983) and smaller localities such as Fântâna lui Mitilan (NISP=139) and La Pietriș (NISP=116) (Supplementary Data 1). There are at least 31 unique taxa from Grăunceanu 11–14; paleoecological reconstructions 15 suggest Grăunceanu was likely an open grassland environment with some woodlands and water sources nearby, comparable to a forest-steppe habitat.

Qualitative analysis of bone surface modifications
We analyzed a total of 4746 specimens from four localities (Grăunceanu = 4524, La Pietriș = 114, Fântâna lui Mitilan = 68, and Fântâna Alortetei = 40). Specimens were examined for bone surface modifications (BSMs) under strong, low-angled light from a gooseneck lamp both macroscopically and with a 10x hand-lens by one of four analysts (SC, BP, SG, CET) following procedures outlined in ref. 17, in constant consultation with each other. Analyses took place during the summers of 2019 and 2022. Fossil surfaces were inspected for a variety of abiotic (e.g., adhering matrix, smoothing, bone surface pitting, erosion/dissolution, cracking/sediment splitting, exfoliation/flaking) and biotic (carnivore or hominin modification, trampling, root/fungal rhizomorph etching, rodent gnawing, insect activity, digestion, preparators marks) modifications and all instances were recorded in a shared Excel spreadsheet. Linear marks, defined as any mark at least twice as long as it is wide, were given particular attention. For any marks we suspected to be cut marks we recorded a subset of the qualitative criteria established in Domínguez-Rodrigo et al. 18, 19: trajectory, presence of barbs, orientation, cross-sectional morphology, number of visible marks, internal symmetry, any shoulder effects visible, any microstriations visible inside the linear mark, any other striations outside of the mark, and color of the interior of the mark relative to the rest of the bone surface. Taxonomic identification, skeletal element, side, portion of element present, and anatomic location of the linear mark (when possible) were recorded. Linear marks of suspected anthropogenic origin were photographed using a DSLR with a macro-lens and a Dino-lite Model Edge and molded two to three times using Coltene President Jet light body dental molding material for further quantitative analysis as described below.

Quantitative analyzes of bone surface modifications. Impressions were taken of marks observed during the qualitative analysis using Coltene President Jet light body dental molding material and sent to MP for quantitative analysis. Marks of interest were molded several times, and the second set of molds were used in this analysis, since adhering sediment on some specimens was removed in the first round of molding. Because all materials are completely mineralized the molding process did not damage or alter the surface of the bones. MP was kept blind to the results of the qualitative analysis and no other contextual information was provided. Photographic images were only shared after the initial analysis when it became clear that some marks analyzed by MP were not the mark intended for study. These were removed from the analysis and in most cases the correct mark was subsequently identified on the mold and analyzed. In some cases, the mark intended for study was not present or was poorly preserved in the mold, and in a few instances marks that were molded correctly were subsequently removed from analysis due to bioerosion overlapping a portion of the mark.

Quantitative analysis of BSMs from Grăunceanu followed the protocol presented in Pante et al. 20. 3D models were created directly from molds using a SNeox 3D optical profilometer (manufactured in 2018) located in the 3D imaging and analysis lab at Colorado State University. All models were produced using a 5x lens that has a z axis resolution of 75 nm. The 5x lens has a numerical aperture of 0.15, a working distance of 23.5 mm, a field of view of 3400 μm x 2837 μm, a spatial sampling of 2.76 μm, and an optical resolution of 0.93 μm. The spatial resolution exceeds that of the Nanovea ST400 white-light non-contact confocal profilometer used in Pante et al. 20 which sampled at 5 μm and 10 μm respectively. However, the z-resolution of 40 nm is slightly higher on the Nanovea profilometer. The scale of these resolution differences does not significantly impact our analysis because they are orders of magnitudes smaller than the scale of the measured differences between mark types. Further, these differences are smaller than the reported variability in measurements taken from a single mark using the same instrument (Pante et al. 20). Processing and analysis of 3D models was carried out using Digital Surf’s Mountains® following Pante et al. 20. Processing included removing outliers, filling in missing data points, and removing the underlying form of the bone with the mark excluded from the form removal process (See Pante et al. 20 for further detail). Data collected through the analysis from the entire 3D model of the BSM were volume, surface area, maximum depth, mean depth, maximum length, and maximum width. Additional data were collected from a profile taken from the deepest point of the BSM including area of the hole, depth of the profile, width, roughness (Ra), opening angle, and radius of the hole.

Fossil data were then statistically compared with a sample of 898 BSMs of known origin, including: 405 cut marks from a variety of stone tool types and raw materials 53; 275 tooth marks from crocodiles and five species of mammalian carnivores 54; 130 trample marks produced by cows on substrates including sand, gravel, and soil 55; and 88 percussion marks from both anvils and hammerstones 56. Surface area and depth of the profile were excluded from the statistical analyses because they are redundant with measurements of volume and maximum depth, respectively. All experimental data were transformed using the Box-Cox method to normalize the distributions for each variable and the same transformations were applied to the archaeological data. Fossil mark classifications were carried out using the quadratic discriminant analysis function from the MASS package 57 in R (version 4.4.1); this method was chosen over a linear discriminant analysis because of unequal covariance matrices and sample sizes in our dataset (following Friedman 58). The accuracy of the QDA model in correctly classifying the experimental BSMs was determined to be 81.4% with 10-fold and 82% with leave-one-out cross-validations performed in R 59 (version 4.4.1) (Supplementary Table 1). Sensitivity values for true positivity, determined using JMP® statistical software (Version 17), range from 0.91–0.99 indicating the model is very successful at classifying all four types of BSMs (see Supplementary Figs. 22, 23, and 24). Prior probabilities were set proportional to the occurrence of each mark type in the dataset (see Supplementary Code 1).

Final determination of cut marks. Any BSM suspected to be anthropic in nature was analyzed both qualitatively and quantitatively. For the qualitative analysis, we recorded the characteristics as described above for linear marks (following 18). Molds of each mark were made and quantitatively analyzed as described above by MP. For these, results included a prior probability that each mark fell into a given category. Some marks were only assessed qualitatively because they were unable to be molded (or scans of the molds failed) or the molds themselves did not pick up the marks in question.

All qualitative and quantitative data were then consolidated and examined collectively by the research team and determinations made regarding the anthropogenic nature of each mark. Linear marks were placed into one of three categories based on their qualitative attributes and results of the quantitative analyses: 1) high-confidence cut marks, 2) probable cut marks, and 3) an “unknown marks” group. High-confidence cut marks were defined as those presenting most of the qualitative characteristics of cut marks as defined by Domínguez-Rodrigo et al. 19 and most were categorized as cut marks in the quantitative analysis. Probable cut marks were those that presented multiple qualitative characteristics indicative of cut marks but were less consistently identified as cut marks in the quantitative analysis or were not analyzed. Unknown marks are those that are not consistent with other taphonomic agents (e.g., trampling, tooth marks) and presented some qualitative characteristics of cut marks but did not include enough characteristics for our team to confidently assign them as anthropogenic in origin. When our qualitative and quantitative analyses returned conflicting results on the nature of the marks, the qualitative analysis was given higher weight due to the holistic nature of our analysis.

U-Pb dating methods. We specifically chose dental samples for U-Pb analysis that would not be useful for other taxonomic assessments (e.g., isolated or fragmentary teeth). Because evidence suggests that the site of Grăunceanu represents a single depositional event, sampling of specimens from this site is essentially random in nature and therefore results should be representative of the entire accumulation. Samples for U-Pb dating were mounted in resin blocks, polished down to reveal internal structures (Supplementary Fig. 25) and arranged in the S155 large-format ablation cell of an Applied Spectra RESOlution-SE ablation system based around an ATL Atlex 193 nm excimer laser. A beam expander was employed for any laser spots above 200 microns. This sample introduction system was coupled to a Nu Instruments Attom-ES high resolution magnetic sector ICP-MS operating in deflector jump mode. Analytical methods follow those outlined for carbonate materials in Woodhead & Petrus 21; details of instrumental parameters are found in Supplementary Table 2.

Laser fluence was typically adjusted to ∼2–3 J cm−2, with a laser repetition rate of 5 Hz. Spot sizes were highly variable in the range 20-200 microns dependent upon sample U contents, which were also highly variable. Tuning of laser gas flows and lenses aimed to achieve maximum Pb counts, without taking 238U into attenuated mode (3 million cps). Low oxide levels (248ThO/232Th < 0.3%) were preferred but no effort was made to tune to Th/U~1 as is common practice in some trace element studies. A brief pre-ablation using a larger spot size was conducted prior to every analysis. Baseline measurement for 30 s was followed by 40 s acquisition during each spot ablation.

Two different primary standards were employed. For most analyses we used a relatively homogeneous ‘in house’ apatite standard ‘BR-13’ but for materials with low U content this was not always possible without tripping the detector into attenuation mode for the calibration standard; for such samples we calibrated using the McClure apatite 60. A variety of apatite secondary standards were employed in the course of this study including Kovdor 61, Durango 62, and ‘401’ 63. All produced ages within 1% of reference values providing confidence in the accuracy of our methods.

Initial data reduction was performed using the VizualAge UcomPbine data reduction scheme (DRS 64) for Iolite, a popular ICPMS data processing software package 65 which is designed to allow the use of heterogeneous reference materials (i.e. those with variable amounts of common Pb such as McClure apatite but assumes no 207Pb/206Pb fractionation). The time-resolved reference material data from a single spot analysis exhibit shifts in 238U/206Pb–207Pb/206Pb isochron space primarily from encountering variably common Pb and/or experiencing downhole Pb/U fractionation. UcomPbine corrects each time slice of background-subtracted data for the reference material based on its known common and radiogenic Pb compositions using a 204Pb-, 207Pb-, or 208Pb-based approach (the 207Pb based correction was used in this study). Drift correction is carried out as usual in Iolite by fitting a function (in this case a smoothing spline) to the reference material analyses that bracket unknowns.

238U/234U measurements were undertaken following the procedure of Hellstrom 66 as modified by Drysdale et al. 67. Initially, relatively large samples of 2–3 mg were drilled from six teeth which returned successful U-Pb isochrons, using a tungsten carbide dental bit. These were dissolved in concentrated HNO3, dried down and re-dissolved in 1.0 ml of 5% HNO3/0.5% HF, 0.01 ml aliquots of which were further diluted in the same acid to 0.8 ml, spiked with a mixed 229Th–233U–236U tracer and equilibrated on a hot plate overnight. These were introduced directly to a Nu Instruments Plasma MC-ICP-MS and measured using a mass jumping routine using mixed Faraday cups and secondary electron multipliers. The procedure was repeated on a second aliquot of each sample, giving inconsistent 230Th/238U ratios but good replication of 234U/238U.

An insoluble organic microparticulate Th-bearing fraction was inferred to be present in the bulk solutions, invalidating 230Th/238U measurements on spiked aliquots of dissolved samples, so all teeth were resampled using a 350 μm tungsten carbide microdrill to a depth of 350 μm giving samples of approximately 75 μg (equivalent to ~10–100 ng U). These were dissolved in concentrated HNO3 and refluxed overnight before addition of 5 drops of conc H2O2 to oxidize the inferred insoluble organic phase, with effervescence observed. Samples were then dried down and redissolved in 0.8 ml of 5% HNO3/0.5% HF, spiked and equilibrated and analyzed as above.

Corrections for initial isotopic disequilibrium in 234U/238U and final isochron construction was performed using DQPB software 68, and the robust isochron fitting technique described by Powell et al. 68 and Monte Carlo uncertainty propagation employing 50,000 iterations.

Stable isotope analysis. Enamel powder from the upper cheek teeth series (P2–M3) of an Equus sp. specimen (VGr.0974) was sampled sequentially along the growth layers using a dental drill, at a resolution of 1.5–2.5 mm, after removing the cementum. Stable isotopes of oxygen and carbon in the structural carbonate of the third molar (M3) were analyzed at the University of Arkansas Stable Isotope Laboratory, USA (UASIL), while the premolars (P2–P4) and the remaining molars (M1–M2) were analyzed at the University of Northumbria, UK. Both laboratories use Thermo Delta V Advantage isotope ratio mass spectrometers coupled to GasBench II sample preparation devices. Prior to analysis, samples were treated with 1M acetic acid, rinsed with deionized water, and dried.

The Northumbria University stable isotope laboratory uses an analysis method adapted from Spötl and Vennemann 69. 2000 ± 20 µg of sample were loaded into 12 mL borosilicate vials capped with butyl rubber septa (LabcoTM), and flushed with helium for 480 s at 1.8 bar on a bespoke flushing box supplied by Sercon Ltd. Carbon. Individual measurement sequences (‘runs’) consisted of 80 measurements, including standards, and all samples were reacted with 102 wt% orthophosphoric acid at a temperature of 70°C, with a reaction time of 104 min prior to measurement to ensure equilibration. For data evaluation, standards were measured at the beginning, end, and approximately every 10 samples in each run. An in-house laboratory calcite standard (Plessen) was used for linearity and drift correction (δ13C=2.40‰, δ18O=−1.31‰, n=10 per run), alongside international standards NBS18 and IAEA 603 for stretching correction (each n≥3 per run). Nominal values from Kim et al. 70 were used for NBS18 (δ13C=-5.01 ‰, δ18O=-23.01 ‰), while values stated by the IAEA were used for IAEA 603 (www.iaea.org; δ13C=2.46‰, δ18O=−2.37‰). One sample of an in-house calcite standard (Pol 2, δ13C=−7.14‰, δ18O=-9.22‰) was added to each run and used to evaluate long-term performance. Based on this secondary standard, for the measurement period relevant to the samples reported here, reproducibility is better than 0.1 ‰ for both δ13C and δ18O. All stable isotope data was initially reported on the VPDB scale, but for climate calculations δ18O values were transferred to the VSMOW scale using the equation δ18O (VSMOW) = δ18O (VPDB) x 1.03092 + 30.9270.

At UASIL, an international standard (NBS19) and an in-house standard (UASIL 22) were used for corrections, each in 8 aliquots analyzed during the entire run. NBS19 standard values are δ13C=1.95‰ VPDB, δ18O=28.72‰ VSMOW, while UASIL 22 values are δ13C=−35.6‰ VPDB, δ18O=13.31 ‰ VSMOW. The standard deviation values for NBS19 and UASIL 22 during the analysis of this set of samples were δ13C–0.09‰ and δ18O–0.17‰, and δ13C–0.08‰ and δ18O–0.13‰, respectively.

Oxygen stable isotopes in the phosphate component of the enamel of M3 were analyzed at UASIL in six samples and reported on the VSMOW scale. The phosphate was precipitated as AgPO4 and analyzed in a Thermo Delta V Advantage IRMS coupled to a TC-EA. Four standards were measured along with this batch of samples: NIST 120c, USGS 80, USGS 81, and an in-house AgPO4 standard, with values of 21.93‰, 13.12‰, 34.78‰, and 6.57‰ VSMOW respectively. Their within-run standard deviation values were 0.60‰, 0.37‰, 0.54‰, and 0.29‰.

Average δ18O values for each tooth were calculated considering that summer season is defined by more data points than winter, thus including all of them in the calculation of the average would skew the result to higher values. To solve this issue, for each tooth we used the same number of samples from summer as there were in winter. For example, P3 has 15 data points from two summer seasons while the full winter season is clearly defined by only 7 points. We eliminated from the calculation of the average the summer data points from both ends of the series until we were left with 7 data points. By contrast, P2 winter is defined by four data points, whereas the summer is defined by only two, leading us to include only two winter values in the calculation.

Using the available δ18O PO4 and δ18O CaCO3 values, we calculated the relationship between them as δ18O PO4 = δ18O CaCO3 x 0.96-8.06 (after rejecting one data point as an outlier), similar to the equation of Iacumin et al. 71 (δ18O PO4 = δ18O CaCO3 x 0.98-8.5). The shift between our measured PO4-CO3 samples is on average 8.9‰, also similar to the Iacumin et al. value of 9.1‰, indicating little to no diagenesis. We therefore have high confidence that diagenetic changes are not meaningfully influencing our results. Because the Iacumin et al. equation is based on a larger dataset, we used that equation to calculate equivalent δ18O PO4 values for the entire data series of δ18O CaCO3 and further, we used the equations δ18O H2O = δ18O PO4 x 0.71+22.60 of Huertas et al. 22, δ18O H2O = δ18O PO4 x 0.9+23 of Kohn and Cerling 24, and δ18O H2O = δ18O PO4 x 1.11-26.44 of Amiot et al. 23 to calculate δ18O values of local meteoric water. Of these equations, only Amiot et al. incorporates error estimates; these errors are given as δ18O w = 1.1128 (±0.0029) x δ18O PO4 - 26.4414 (±0.0508). However, integrating those uncertainties would affect only the second digit of the results, with no discernible change in our interpretation.

While enamel δ13C values are generally shifted by 14.1±0.5‰ compared to plants 72, we calculated the average annual plant value and corrected it for the difference of 1.7‰ between modern and Quaternary atmospheric CO2 73. The resulting plant tissue value was then used in the equation of Kohn 29, which allows for the calculation of the mean annual precipitation amount (MAP): δ13C (‰, VPDB) = −10.29 + 1.90 × 10-4 × altitude (m) − 5.61 × log (MAP+300, mm/year) - 0.0124 x latitude (°), where the altitude of our site is 300 m and latitude is 45°. Kohn and McKay 74 estimate the uncertainties of this equation as being ~50% for MAP values higher than 500 mm/year, but never decreased below 120 mm/year. Although Kohn 29’s equation is intended to calculate annual precipitation values, we divide the results into 12 months in order to calculate monthly values (as described in the SI). To our knowledge, such use of this equation has not been validated before.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Curran et al discuss evidence suggesting the presence of hominins in Romania approximately 2 million years ago, predating the widely accepted oldest Eurasian site of Dmanisi (Georgia). Their evidence includes cut-marked bones and new Uranium-Lead (U-Pb) dates of fossil dentine. The authors also use horse enamel isotope data for a coarse paleoclimate and paleoenvironmental reconstruction. Except for a few typos and awkward sentences, the manuscript is generally well- written. Figures need a little extra work (see below), but they are fairly well-crafted. Unfortunately, I must defer the evaluation of the hominin presence and U-Pb components of the paper to other reviewers, as they are outside my area of expertise. I can only say that their arguments and conclusions seem reasonable in light of the data to a non-expert.\n\nMy main concern with this manuscript is the oversimplification of the interpretation of the isotope data. More specifically, here are the two main issues:\n1) Interpretation of enamel oxygen isotope data. The authors interpret lower-than-present precipitation d18O (calculated from horse enamel d18O) as an increase in winter precipitation component to aquifer recharge relative to the present. This is one of the many possible explanations. Alternative possibilities include a lower mean annual temperature, and a change in air mass trajectories. These alternative possibilities must be discussed here in depth. Furthermore, there is no evidence that the horses were drinking groundwater. Finally, the uncertainty in the average d18O value of precipitation is largely underestimated. The authors need to perform a rigorous error analysis and propagate the errors from ALL sources.\n2) Interpretation of enamel carbon isotope data. There are many issues with the calculation of MAP from enamel d13C, including a) the possibility that the investigated taxon was selective in terms of its diet and habitat; b) the uncertainty in the enrichment factor between diet and enamel; and c) the uncertainty in the relationship between the d13C of modern C3 plants and MAP. Like with d18O, the authors largely underestimate the error in the reconstructed MAP. I am inclined to say that because of all these issues these MAP calculations and paleoenvironmental reconstructions have such a large uncertainty that they are almost meaningless, especially if only one taxon is used like in this case. The authors should consider eliminating the calculated MAP from the manuscript unless they provide a robust estimate for the uncertainty.\n\nThese two issues must be fixed before the manuscript can be considered for publication in Nature Communications. Minor issues include hard to read figures and text. For instance, the text in Fig. 1 is hard to read (especially in the inset) and the “18” and “13” in Fig. 2 should be superscripts. Again, the manuscript also needs some minor editorial work to fix some typos and to rewrite some poorly structured sentences.", "reviewer_2": "The paper by Curran and colleagues titled “Hominin presence in Eurasia by at least 1.95 million years ago” presents a highly important find that hominins were indeed present and active outside of Africa at least 1.95 million years ago. I have thoroughly enjoyed reviewing this paper, and certainly think it is relevant and falls within the scope of this journal. As the authors themselves point out, however, extraordinary claims require extraordinary evidence to avoid and answer possible scepticism. While I think the photographic evidence of the cut marks is very compelling, I will focus on assuring that the classifications and analyses are as rigorous as possible. For this reason, I will therefore dedicate my review solely to the taphonomic questions of the paper, as all other aspects are far from my area of expertise.\n\nFirst of all, I think a terminological clarification is needed, as I find the use of the word score throughout the text to be quite problematic. The word “score” typically refers to carnivoran activity, as in “tooth score”, while I have very rarely seen people refer to cut marks as scores. It is fundamental that we are as clear as possible in the words we use to describe certain processes, so as to avoid confusion, and because of this I encourage the authors to use a different term that is less likely to be mistaken with other taphonomic alterations (we already have enough issues with taphonomic equifinality, we don’t want to add any “terminological equifinality” to the mix!). In the taphonomic atlas published by Ferndández-Jalvo and Andrews, for example, the authors talk about a “linear” category, where both tooth scores and cut marks can be considered “linear marks” – they have a length that is more than twice their width. I think it is much more appropriate to use a term like this, than use a word that is already shared in the name of a different type of taphonomic alteration.\n\nI also think more evidence about how the collections were obtained and how they are currently being curated is also very important, which ties in to my next comment about the moulding process. The authors state that at least two moulds were taken per cut mark, however a recent study by Valtierra et al. (2023, DOI: 10.1111/arcm.12938) have shown that successive extraction of moulds can have an effect on the morphology of the traces being studied. While they showed that statistically these changes are not too important, this work was carried out on the morphological study of a single profile along the BSM, at a very specific point. It is also important to point out that with each mould taken, the natural variability of metric variables decrease, which could have a potential impact on certain analyses. While some of the variables used in the quantitative analysis of the present study are similar and should not, therefore, change too much, I am pretty sure that there will be a change in the morphology of the trace after each mould is taken. This leads me to ask, which moulds were used for the quantitative analysis? The first mould taken? Or the last mould taken? How were the fossil specimens cleaned as well? When performing such high resolution analyses, and considering variables such as volume, these different processes could have an impact. While I don’t think this would have too much of an impact, it is important to at least document and be aware of it.\n\nBeyond overlying trephic modifications, I also think it is very important to consider the impact dendritic etchings and root marks may have had on the samples, especially when considering the methodology employed for quantitative BSM analysis and identification. For example, in VGr.0519, VGr.1557, VGr.1767, VGr.2275, VGr.2687 and VGr.2649 you can see overlapping biochemical alterations that would surely affect the calculation of BSM volume and other parameters, wouldn’t these parameters now be unreliable for the identification of the traces? To what extent should the analysis therefore include these variables? What happens when you perform calculations on some of the samples including or excluding some of these parameters? Do the classifications change? It is very true that almost any methodology dedicated to the high- resolution analysis of BSMs would be conditioned by these alterations, however, calculations such as volume are more likely to be effected by such changes to the cortical surface.\n\nWith regards to data acquisition, much more information is needed regarding the microscopic analysis, including the precise FOV and numerical aperture of the lens used and the windows scanned; variables that have been proven in the past by authors such as Calandra et al. (2019, DOI: 10.1038/s41598-019-42713-w) to effect the replicability of analyses. Also, if I recall correctly, the S Neox has the option of obtaining 3D models using either Focus Variation or Confocal techniques, that would also condition the quality of the 3D model produced, so the modality must be stated. Why was only 5x magnification used? What colour light was used? My own experience with a similar microscope found 5x and 20x to produce different results for parameters such as roughness and volume when doing use-wear analysis. For this reason I’d also argue that comparing data obtained from the S Neox and the Nanovea ST400 could definitely be considered problematic for some of the quantitative variables included in the study (again, see the work of Ivan Calandra). The authors state that the results should not be impacted because of “coarser scales of measurement” (which would probably answer my question of why was only 5x used), however this is yet to be tested, and when you calculate variables such as volume and roughness, then multiple researchers, including myself, would in fact disagree with this statement. Information about how many non-measured points during the scans are also required information, so as to properly assess the quality of the scans used for the analysis.\n\nRegarding the post-processing of the 3D models obtained using the microscope, the author states that they removed form from the data, yet no indication about what polynomial was used to do this. I have to say that I think form removal for cut mark analysis is potentially very problematic, however I understand why it has been done considering some of the variables that have been measured. Form is removed in microwear studies because this directly has an impact on roughness parameters, and will also influence how easy it is to compute each traces volume, however when we are talking about the properties of the profile, form removal will have an effect on mark morphology.\n\nWith regards to the statistical analysis of the data, and classification tasks performed, I also have some concerns that need addressing. First of all, Quadratic Discriminant Analysis (QDA) should not be used for ordination, as has been done for supplementary figure 22. While use of QDA algorithms is perfectly acceptable for classification, it should not be used for ordination for a number of mathematical and epistemological reasons. QDA is a classification algorithm, so of course there will be separation between samples in the constructed and visualised feature spaces, because that’s what it is designed to do – force separation to exist. There is also a huge issue with these types of constructed feature spaces due to issues with sample size and imbalance. While the individual to variable ratio is acceptable for the dataset provided, the sample imbalance is large enough to be problematic, in particular between cut marks and percussion marks (side note – why are these two compared in general, are there people who consider equifinality to exist between percussion and cut marks?). What the authors should use instead, if they wish to visualise the multivariate differences between samples, is perform and visualise a PCA with a biplot, or simply present boxplots for each variable separately. If they want to use the QDA figure, then it must be used as a complementary representation of the data. Its also important to ensure that the feature space is constructed using cross-validation, otherwise separation is much more likely to be spurious in nature. For an example, see https://doi.org/10.1007/s11692-020-09494-x - although this paper was done for the purpose of discussing use of similar techniques in geometric morphometrics, the algebraic properties of the ordination technique are very similar, while the input to the figures in this paper are not landmark data but instead standard multivariate distributions and are thus completely comparable to the present case study. Finally, an ordination technique like these ones should not be presented without including in the axes labels the percentage of informative that is captured by each dimension, otherwise they are very hard to interpret.\n\nFor the evaluation of how well a classification algorithm works before being applied to archaeological data (presented in the first paragraph of page 22), the “accuracy” metric is often misleading, especially in cases where the imbalance is notable. Instead the authors must present other metrics such as precision, recall, sensitivity or specificity, that are more informative of other parts of the confusion matrix\n\nSome minor comments:\nIn the results section of the main article the authors state that a reduced number of bones presenting polishing is proof of limited transport. Polish does not necessarily imply movement, it’s the action of water-born abrasive particles altering the overall surface of the bone, this could occur with a low flowing course of water over a prolonged period of time without moving the bones at all. In the supplementary materials the authors mention rounding as well. I think it’s more informative to simply present all of the water-based alterations to cortical surfaces in the main text as well.\n\nIn the last paragraph of the taphonomy section in the results the authors mention that some specimens have “more degraded surfaces”, it would be nice if the nature of this degradation be clarified – is it due to fluvial processes, for example?\nHave percussion marks or evidence of anthropogenic fracture patterns been identified in the site?\n\nOn line 227 the authors include some citations of papers that support the notion that sites exist with identified cut marks but no associated lithic assemblages. Among these cited papers, they include BK from Olduvai, however BK does contain a huge amount of lithics published as far back as the Leakeys, and is even mentioned in the abstract of the very paper they cite. I can assure you BK has lithics, while Mary Leakey even noted the presence of handaxes from this site. You also cite Dikika, which is very problematic, and many would even argue that not only are there no lithics, but no cut marks either\n\nIn summary, I once again wish to thank the authors for their work, I truly have no problem at all with the validity of their finds, I just hope to ensure that the presentation of results are as rigorous, robust, and transparent, as possible. I have no problem disclosing my name as a reviewer - Lloyd A. Courtenay", "reviewer_3": "The paper presented by Curran and colleagues, entitled \"Hominin presence in Eurasia by at least 1.95 million years ago,\" presents clear evidence of ancient anthropogenic activity in the form of butchery marks from several sites in the Olteț River Valley in Romania. The remains come mainly from a single Grăunceanu assemblage, although they also provide data from some other penecontemporaneous sites. A list of U-Pb dates establishes an age of 1.95 Ma. This work offers highly relevant information for understanding human dispersal prior to 2 Ma and provides palaeoenvironmental data, through isotopic analysis, suggesting conditions conducive to the development of such human movements.\n\nThe authors present bone remains with unquestionable cut marks, which have been analyzed using both quantitative and qualitative methods. Nevertheless, it is surprising that the authors state that the fossils were examined under strong, low-angled light for evidence of bone surface modifications. It would have been more appropriate to use microscopic observational methods (stereomicroscope, ESEM, digital microscopes, etc.). It seems plausible to suggest that this is the reason why the analysis method of Rodrigo et al., 2009, has been modified, with the aim of avoiding the microscopic variables of the analysis. It is likely that had a microscopic approach been employed, a greater number of marks would have been identified. The marks attributed to cut marks are clearly visible to the naked eye, and therefore it is reasonable to hypothesize that a microscopic analysis would have identified a greater number of evidence, as well as helped to attribute some of the probable marks to high-confidence marks.\n\nThe authors propose that some of the striae are attributable to high-confidence cut marks, while others may be associated with probable or unknown marks. Although they employ a quantitative methodology, they utilize a similarity-by-association approach to ascertain the origin of the marks as either anthropogenic or non-anthropogenic. It would be preferable for the authors to not use association to attribute the marks, but rather to only apply the proposed analytical principles.\n\nIt is unclear why the teeth were excluded from the taphonomical analysis. If one of the objectives was to identify butchering marks, the teeth may show evidence of cuts related to the cheek defleshing and/or tongue removal. In addition, some of them could have been dated, obtaining direct dates from remains modified by lithic tools.\n\nIn the text (third paragraph of the Taphonomy section), the number of remains with high-confidence cut marks is initially indicated as seven. However, later in the text, it is stated that there are eight remains, and these are also the number of fossils that appear in Table S2. Therefore, it is necessary to correct the text.\n\nIn some instances, the observed marks appear to be clearly anthropic in origin. However, the analysis suggests otherwise. This discrepancy may be a consequence of the macroscopic analysis developed by the authors and the omission of certain microscopic variables that could have biased the results. A good example is the metopodial VGr.0429. Upon examination of the image (Fig. S12), it appears to be a clear anthropic disarticulation mark. However, the analysis suggests that it is, in fact, a tooth mark. It is recommended that these probable cutmarks be reanalyzed under microscopic tools (a stereomicroscope, at a minimum).\n\nVGr.1490 is described as a cutmark with a high probability and a high confidence, yet it appears both in the text (page 6 of SI) and in table S2 as a probable cut mark. Could this be an error? In the SI description, the authors even assign the activity carried out by hominins, tendon extraction.\n\nWhy do the authors attribute the mark observed on specimen VGr.2052 (Fig. S18) to trampling if the analysis indicates that it is a cut mark with a high probability? This suggests that the authors sometimes adhere to the result of the quantitative analysis, but at other times, the observation of the researchers takes precedence.\n\nGiven the paucity of information on the stratigraphic context of these ancient collections, it is understandable that the authors refrain from offering interpretations regarding the type of occupation or the mode of access to the resources. However, it is crucial to clarify the authorship of the marks and their age.\n\nTables 1 and S7 are invaluable for the compilation of all ancient sites with evidence of hominids or anthropogenic activity (lithics and butchery marks) and the corresponding bibliography. It is recommended to include additional bibliographical references concerning published extended taphonomic studies of the cited sites:\n- For the Ain Boucherit site (Algeria):\nCáceres, I., Chelli Cheheb, R., van der Made, J., Harichane, Z., Boulaghraief, K., Sahnouni, M., 2023. Assessing the subsistence strategies of the earliest North African inhabitants: evidence from the Early Pleistocene site of Ain Boucherit (Algeria). Archaeological and Anthropological Sciences 15, 87. https://doi.org/10.1007/s12520-023-01783-8\nFor the Sima del Elefante site (Spain):\nBlasco, R., Blain, H.A., Rosell, J., Díez, J.C., Huguet, R., Rodríguez, J., Arsuaga, J.L., Bermúdez De Castro, J.M., Carbonell, E., 2011. Earliest evidence for human consumption of tortoises in the European Early Pleistocene from Sima del Elefante, Sierra de Atapuerca, Spain. Journal of Human Evolution 61, 503–509. https://doi.org/10.1016/j.jhevol.2011.06.002\nHuguet, R., Saladié, P., Cáceres, I., Díez, C., Rosell, J., Bennàsar, M., Blasco, R., Esteban-Nadal, M., Gabucio, M.J., Rodríguez-Hidalgo, A., Carbonell, E., 2013. Successful subsistence strategies of the first humans in south-western Europe. Quaternary International 295, 168–182. https://doi.org/10.1016/j.quaint.2012.11.015.\n\nThe evidence presented is of significant interest and should be published. The utilization of microscopy to analyze the attributes of the cut marks is recommended in order to identify additional marks, correctly assign probable ones, and enhance the images presented. It is also recommended to add some high-resolution images of examples of trampling or carnivore toothmarks.", "reviewer_4": "Hominin presence in Eurasia by at least 1.95 million years ago, by Curran et al.\nThis is a very interesting work with some exciting results suggesting early evidence of human presence in East Europe by > 2 Ma, i.e. much older than thought until very recently. Please note that the evaluation of the cut marks on bones and whether there are human made is beyond my area of expertise. I will mostly focus my review on the chronology and general context of the fossils.\n\nAs a general comment, I would say that the manuscript is overall well written and easy to follow. The results and data presented apparently support the interpretation proposed by the authors. Actually, the main issues are in my opinion not related to what is shown in the work, but rather to what is NOT shown. In that regard, the absence of a detailed contextualization of the fossils is my major concern. I get that the fossils are coming from old excavations and that the various sites are most likely no longer accessible, but we have virtually no idea about their stratigraphic position, and whether the fossils represent one single or various assemblages. In other words, the integrity of the fossil assemblages is simply unknown. This major question mark about the integrity of the fossil assemblage, and whether all fossils correspond to a single depositional event, is actually not addressed in the present work, since the fossils dated with U-Pb are not those that have cut marks. While one may assume that they have all the same age, we actually don’t know, and this should be verified by running additional analyses on the cut-marked bones. The same comment applies to the stable isotope analyses performed on the equid teeth. Direct U-series dating should be performed in order to confirm its antiquity and make sure it can be closely associated to the other fossils. In that regard, I genuinely wonder whether Rare Earth Elements could also help demonstrate that all fossils have the same origin and come from the same stratigraphic unit.\n\nThe whole context is briefly summarized in Methods, but remains (too) vague. The authors mention that the location of various localities have been ‘reidentified’, therefore it should be possible to get further information about the context, and establish a proper stratigraphic column, description of the sediment, with some additional information about the formation of the fossil horizon and the depositional environment for the 3 localities. A geological or geomorphological map of the ORV with the stratigraphic correlations of the various sections would also be very useful.\n\nThe lack of independent age control is also another major concern. There is some circular reasoning in using the biochronology (which is based on independently dated assemblages) to support the numerical ages obtained here. The uncertainty on biochronological inferences is usually minored and not really taken into account, and especially on large mammals. I doubt very much that the precision on the biochronological age is sufficient to undoubtedly position the site around 2.2-1.9 Ma, within a small 0.3 Ma age window. This would need further explanation with a dedicated section in SM, and in particular about the robustness of the age control on the localities used as chronological tie points. Rodents may possibly help refining the chronology of the assemblage (e.g. Maul & Markova QI, 2007), but seems like from Table S1 that there are too few remains. As mentioned above, the absence of independent age control adds a non-negligible uncertainty to the dating results, especially when knowing that dental tissues behave as open system, and may show U leaching, which would artificially age the sample. I had a quick look at Andreescu et al. QI (2013), and I didn’t find the existing age control on the Tetoiu Fm very convincing to securely confirm its chronological attribution to the Early Pleistocene.\n\nOn an unrelated matter, I see in Table S5 that several samples are in secular equilibrium, while other do show U-leaching (i.e., beyond equilibrium values). For example, this is the case for sample VGr.2256-2. My first reaction would be that this outcome may seriously question the reliability of the U-Pb age of 1.95 Ma obtained on that same tooth. I think this aspect would deserve some additional discussion in the manuscript. From the information given in SM, I understand that both solution and laser ablation analyses were carried out on the teeth. Some samples were microdrilled. I am not sure to understand why, and this could be clarified in the text perhaps. It is also unclear to me how LA analyses were performed, either by drilling vertically within the dentine, or through mapping an area of the occlusal surface? If so, can we see the maps with the spatial distribution of the various ratios? Please add some extra explanations in SM. The exact location of the samples analyzed for U-series (given in Table S5) within the teeth could also be given. Fossil teeth are known to be highly heterogeneous in terms of spatial distribution of radioelement, and especially old teeth (e.g., Duval et al. QG, 2011) and it is unclear to me how this uncertainty was evaluated. Some domains within the dentine or enamel may locally show very distinct U-series data, which would bias the age results. It would be strongly advisable to sample other parts of each tooth to evaluate this variability and help identify more recent uranium uptake or leaching overprints.\n\nSpecific comments\n1. Main text, ‘The current best’: What does ‘best’ mean here? Please be more specific.\n2. What is the definition of ‘ephemeral traces’ here? There several occurrences of this term throughout the manuscript and I guess this should be explained somewhere.\n3. U-Pb dating, ‘2.01±0.20 Ma to 1.87±0.16’. Please indicate at least once in the main text whether errors given are 1 or 2σ. ‘Ma’ missing after ‘1.87±0.16’\n4. U-Pb dating, ‘which is ubiquitous in fossil teeth and bone samples’: Do you mean at this specific site? Please clarify the sentence. There are actually many works showing the timing of the U uptake may significantly vary within a given dental tissue or bone (e.g., Grün and Taylor Ancient TL 1996; Grün, RM 2009; Duval et al., QG, 2011), between the various tissues of a given tooth, and between fossils from the same site/horizon.\n5. U-Pb dating, ‘Grăunceanu took place >1.95 Ma and potentially >2.01±0.20 Ma’: actually these two ages mentioned here are virtually the same given the associated error. There is no need in my opinion to mention them both, one is sufficient.\n6. U-Pb dating, ‘These results make these sites the most securely dated localities from Early Pleistocene Europe’: the concept of ‘securely dated’ is very relative. The absence of independent age control and the presence of minimum ages on fossils do not make the sites better constrained than others in my view.\n7. Discussion, ‘Our paleoenvironmental data for Grăunceanu indicate a more open, arid environment with some nearby woodlands and water sources’ : do you mean more open and arid than today? Please specify. Additionally, the estimated MAP for Grăunceanu is higher than current conditions, which seems to be in contradiction with the statement above on a more arid environment.\n8. Discussion, ‘While it is clear that the hominin presence in Eurasia at this time was not extensive’: Why is it clear? And what do you mean by extensive? There is actually early evidence of human presence in Europe from Ukraine to Spain.\n9. Fig. 3 caption: please specify whether errors are displayed at 1 or 2σ .\n10. Fig. S24: I am not sure to understand the interest of showing elemental concentrations in dental tissues? What kind of information does it bring?\n11. Fig. S24: Metal concentration or elemental concentration?\n12. Please would it be possible to show instead the maps of U concentrations, U234/U-238 and apparent age for each of the dated samples?\n13. Fog. S25: How was the present-day U234/U-238 measured?\n14. Table S1: Does the composition of the fossil assemblage come from previous studies (if so, please cite the relevant references), or from the direct study of the specimens by the authors? If the data given for this table are coming from a new taxonomic revision of the fossil assemblage, then a dedicated section including a proper systematic paleontology study should be included.\n15. Table S4: which U-234/U-238 values were used for the disequilibrium corrected age? This may be clarified, as Table S5 shows that there are several values available for each sample.\n16. Table S5:\n- Please clarify whether MXF and UMF samples have all been measured through solution analyses, or UMF samples have been analyzed using Laser Ablation?\n- Where (=which lab) was Th-230/Th-232 ratio measured?\n- Not sure I understand well the U mobility column. Do ‘Closed’ mean ‘secular equilibrium’, and ‘U gain’ mean that the ratio are still below secular equilibrium?\n- The apparent U-Th age can be added to this table. Out of curiosity, how ca we explain that the same tooth may provide an apparent U-Th age of <500 ka and a U-Pb age of > 1.9 Ma. This would probably deserve some additional comments in SM.\n17. I have my doubts on the reliability of the age of many localities mentioned in that table. Can the authors at least add one column mentioning the dating methods employed to constrain the age of the site (e.g., U-Pb, palaeomagnetism, biochronology, ESR, TCN, etc.)?"}, "author_response": "We greatly appreciate the time and effort the four reviewers put into their suggestions, which\nhave improved this manuscript. In order to make a more readable revision, we have made a\nnumber of changes to align with Nature Communications specific format that were not tracked\n(formatting of headings, supplementary figure references, etc.). All changes to the manuscript\n(and supplemental documents) to address reviewer comments have been tracked.\nBelow we provide a point-by-point (bolded) response to the reviewers’ (non-bolded) comments.\nPlease do not hesitate to reach out with further questions and clarifications.\nReviewer #1 (Remarks to the Author):\nCurran et al discuss evidence suggesting the presence of hominins in Romania approximately 2\nmillion years ago, predating the widely accepted oldest Eurasian site of Dmanisi (Georgia). Their\nevidence includes cut-marked bones and new Uranium-Lead (U-Pb) dates of fossil dentine. The\nauthors also use horse enamel isotope data for a coarse paleoclimate and paleoenvironmental\nreconstruction. Except for a few typos and awkward sentences, the manuscript is generally well-\nwritten. Figures need a little extra work (see below), but they are fairly well-crafted. \nUnfortunately, I must defer the evaluation of the hominin presence and U-Pb components of the\npaper to other reviewers, as they are outside my area of expertise. I can only say that their\narguments and conclusions seem reasonable in light of the data to a non-expert.\n\nMy main concern with this manuscript is the oversimplification of the interpretation of the\nisotope data. More specifically, here are the two main issues:\n1) Interpretation of enamel oxygen isotope data. The authors interpret lower-than-present\nprecipitation d18O (calculated from horse enamel d18O) as an increase in winter precipitation\ncomponent to aquifer recharge relative to the present. This is one of the many possible\nexplanations. Alternative possibilities include a lower mean annual temperature, and a change in\nair mass trajectories. These alternative possibilities must be discussed here in depth.\n--We have now clarified how we came to this conclusion by adding estimates of winter and\nsummer precipitation amounts, which we summarize in the main text and then detail in the\nsupplementary discussion. Regarding potential changes in air mass trajectory compared to\nthe present, we do not have information about the Early Pleistocene setting in this region.\nAt present, two main, isotopically different vapor sources are at play in Romania, namely\nthe North Atlantic and the Eastern Mediterranean (which contributes water vapor with\nhigher isotopic values due to increased evaporation at source). Local meteoric water lines in\nthis region have equations similar to the global meteoric water line, indicating the North\nAtlantic as being the most important moisture source (Dragusin et al., 2017). We have\nadded some text that mentions these possible influences but note that we either consider\nthose possible interpretations unlikely due to the fauna at the site (lower mean annual\ntemperatures), or unable to incorporate these sources into our analyses because of\nuncertainties regarding the fossil record (changes in air mass trajectories).\nFurthermore, there is no evidence that the horses were drinking groundwater.\n1\n--We do not suggest that horses drank groundwater. In the Supplementary materials, we\nwrite: “Horses, which are obligate drinkers (water dependent), take most of their water\nfrom springs, rivers, or lakes, each of these sources being defined by their own isotopic\ndrivers such as seasonality of recharge, evaporation, mixing, or recharge altitude. Rivers\nhave lower amplitude of isotopic variability than rainfall, as they drain aquifers that are\nmost of the time well mixed throughout the year.” We have clarified this in the main text.\nFinally, the uncertainty in the average d18O value of precipitation is largely underestimated. The\nauthors need to perform a rigorous error analysis and propagate the errors from ALL sources.\n--We recognize that there are four potential sources of error when reconstructing δ18O of\nmeteoric water: 1) analytical measurement of samples, 2) the PO4-CO3 relationship, 3)\ndiagenesis, and 4) meteoric water to enamel transfer functions. We have tried to clarify\nthese complicating factors in the text and we address them here as well.\n• Analytical errors are stated as 0.1‰ for CO3 at both labs, while it was 0.60‰,\n0.37‰, 0.54‰, and 0.29‰ for PO4, for the four standards used.\n• We directly assessed the PO4-CO3 relationship for our M3 sample; this relationship\nis defined as δ18OPO4 = δ18OCaCO3 x 0.96- 8.06, which is very similar to that of\nIacumin et al. (1996), δ18OPO4 = δ18OCaCO3 x 0.98 - 8.5. The shift between our\nmeasured PO4-CO3 samples is on average 8.9‰, also similar to the Iacumin et al.\n(1996) value of 9.1‰, indicating little to no diagenesis. We therefore have high\nconfidence that diagenetic changes are not meaningfully influencing our results.\n• Though we consider our original analyses to be robust, based on the reviewer’s\nconcerns regarding sources of error we recalculated our δ18OPO4 estimates using\nthe equation of Iacumin et al. (1996) instead of our own equation. Although the\nIacumin et al. (1996) PO4-CO3 relationship is not provided with uncertainties, we\ndecided to use it because it is based on a larger dataset, instead of using our own\nrelationship derived from only five sample points. In addition, their equation also is\nlikely to be more robust because they have smaller measurement uncertainties (0.15-\n0.20‰). In contrast our measurement run had higher values, which would render\nany error propagation from our equation meaningless. However, there are no errors\nto propagate from the Iacumin et al. (1996).\n• Of the three equations we then employed to reconstruct meteoric isotopic values\n(Delgado Huertas et al., 1995; Kohn and Cerlin, 2002; Amiot et al., 2002), only\nAmiot et al. incorporates error estimates. These errors associated with δ18O in\nenamel PO4 and meteoric water are given as δ18Ow= 1.1128 (+/- 0.0029) x δ18OPO4 -\n26.4414 (+/- 0.0508), and integrating those uncertainties would affect only the\nsecond digit of the results, with no discernible change in our interpretation.\nFurthermore, all three equations (Delgado Huertas et al., 1995; Kohn and Cerling,\n2002; Amiot et al., 2002) return mostly similar results and by averaging their results\nfor the P2-M3 series we obtained the value of -10.72 (+/-1.02‰), which is still below\npresent-day values.\n2\n• To conclude, our PO4-CO3 relationship resulting from a small dataset is similar to\npublished values based on a much larger dataset of Iacumin et al. (1996) and shows\nno important signs of diagenesis. Out of three published transfer functions from\nPO4 to meteoric water only one is provided with uncertainties and those are very\nlow, and all three return similar values, resulting in an average value with a\nstandard deviation of 1.02‰. Thus, a full propagation of errors as outlined in Amiot\net al. (2002) only has a minor effect on the final values. We therefore consider the\nstandard deviation of 1.02‰ to be the most important statistical parameter for\ninterpreting our results; we have added this standard deviation to Figure 4 to make\nthis clearer. We have also edited the text to reflect the change in our use of Iacumin\net al.’s equation.\n2) Interpretation of enamel carbon isotope data. There are many issues with the calculation of\nMAP from enamel d13C, including a) the possibility that the investigated taxon was selective in\nterms of its diet and habitat;\n--The reviewer is correct that our investigated taxon, Equus, is highly selective in its diet\nand habitat, but that is precisely why we feel this is a powerful taxon to investigate. Modern\nEquus preferentially graze throughout the year, and all existing data suggest this has been\nthe case for this genus since the Pleistocene. Because this taxon is so restricted in its diet,\nwe have especially high confidence that we are actually capturing seasonal changes in mean\nannual precipitation (and/or migratory behaviors), rather than seasonal dietary shifts.\nOther analyses by our research team (Curran et al., 2021) have investigated isotopic\ndietary signals in ruminant taxa (i.e., deer), as well as meso- and microwear analyses of the\nteeth to get at dietary behaviors in species from Graunceanu. As we state in the\nsupplementary materials discussion of the stable isotopes: “Plant tissue δ13C reconstructed\nfrom enamel δ13C values varies between -28.9‰ ±0.5‰ and -25.7‰ ±0.5‰ (VPDB), which\nis consistent with feeding in a woodland to woodland-mesic C3 grassland (as defined by\nDomingo (2013), and supports prior environmental reconstructions for Grăunceanu\n(Curran et al., 2021).”\nb) the uncertainty in the enrichment factor between diet and enamel;\n--The isotopic enrichment between diet and enamel is given by Cerling and Harris (1999) as\n14.1 +/- 0.5‰ and we integrated the uncertainty in the new calculation of the seasonal\nprecipitation amounts. This results in precipitation amount estimates for the driest month\nof 6-22 mm/month, while for the wettest month they are between 90 and 148 mm/month.\nThe text has been updated to include these new ranges of precipitation.\nand c) the uncertainty in the relationship between the d13C of modern C3 plants and MAP. Like\nwith d18O, the authors largely underestimate the error in the reconstructed MAP. I am inclined to\nsay that because of all these issues these MAP calculations and paleoenvironmental\nreconstructions have such a large uncertainty that they are almost meaningless, especially if only\none taxon is used like in this case. The authors should consider eliminating the calculated MAP\nfrom the manuscript unless they provide a robust estimate for the uncertainty.\n--We appreciate the reviewer’s concerns regarding our MAP calculations and their\nuncertainties. We agree that our previous estimates of MAP had a broad range (868-1134\nmm/year) even though they did not include the 0.5‰ uncertainty from Cerling and Harris\n3\n(1999); including those uncertainties would certainly increase those ranges in such a way\nthat they are meaningless. Unfortunately, the uncertainty in the relationship between the\nδ13C of modern C3 plants and MAP as given by the equation by Kohn (2010) is unknown\nsince they did not provide uncertainties. Thus, instead of taking our original approach, we\nnow instead discuss the reconstructions for just the driest (i.e., summer) and wettest (i.e.,\nwinter) months in order to define seasonal variability in precipitation. The text has been\nedited to reflect this change in focus.\nThese two issues must be fixed before the manuscript can be considered for publication in\nNature Communications. Minor issues include hard to read figures and text. For instance, the\ntext in Fig. 1 is hard to read (especially in the inset) and the “18” and “13” in Fig. 2 should be\nsuperscripts. Again, the manuscript also needs some minor editorial work to fix some typos\nand to rewrite some poorly structured sentences.\n--Figure 1 will be submitted as a higher resolution image and Figure 4 (the isotope figure)\nhave been edited.\n--We also have gone back through the manuscript to correct any typos we found and\nchecked for awkward sentences, though we cannot be certain which ones Reviewer 1 is\nreferring to specifically.\nReviewer #2 (Remarks to the Author):\nThe paper by Curran and colleagues titled “Hominin presence in Eurasia by at least 1.95 million\nyears ago” presents a highly important find that hominins were indeed present and active outside\nof Africa at least 1.95 million years ago. I have thoroughly enjoyed reviewing this paper, and\ncertainly think it is relevant and falls within the scope of this journal. As the authors themselves\npoint out, however, extraordinary claims require extraordinary evidence to avoid and answer\npossible scepticism. While I think the photographic evidence of the cut marks is very\ncompelling, I will focus on assuring that the classifications and analyses are as rigorous as\npossible. For this reason, I will therefore dedicate my review solely to the taphonomic questions\nof the paper, as all other aspects are far from my area of expertise.\nFirst of all, I think a terminological clarification is needed, as I find the use of the word score\nthroughout the text to be quite problematic. The word “score” typically refers to carnivoran\nactivity, as in “tooth score”, while I have very rarely seen people refer to cut marks as scores. It\nis fundamental that we are as clear as possible in the words we use to describe certain processes,\nso as to avoid confusion, and because of this I encourage the authors to use a different term that\nis less likely to be mistaken with other taphonomic alterations (we already have enough issues with\ntaphonomic equifinality, we don’t want to add any “terminological equifinality” to the\nmix!). In the taphonomic atlas published by Ferndández-Jalvo and Andrews, for example, the\nauthors talk about a “linear” category, where both tooth scores and cut marks can be considered\n“linear marks” – they have a length that is more than twice their width. I think it is much more\nappropriate to use a term like this, than use a word that is already shared in the name of a\ndifferent type of taphonomic alteration.\n--We have changed the generalized form of “scores” to “linear marks” throughout the\nmanuscript.\nbeing curated is also very important, which ties in to my next comment about the moulding\nprocess. The authors state that at least two moulds were taken per cut mark, however a recent\nstudy by Valtierra et al. (2023, DOI: 10.1111/arcm.12938) have shown that successive extraction\nof moulds can have an effect on the morphology of the traces being studied. While they showed\nthat statistically these changes are not too important, this work was carried out on the\nmorphological study of a single profile along the BSM, at a very specific point. It is also\nimportant to point out that with each mould taken, the natural variability of metric variables\ndecrease, which could have a potential impact on certain analyses. While some of the variables\nused in the quantitative analysis of the present study are similar and should not, therefore, change\ntoo much, I am pretty sure that there will be a change in the morphology of the trace after each\nmould is taken. This leads me to ask, which moulds were used for the quantitative analysis? The\nfirst mould taken? Or the last mould taken? How were the fossil specimens cleaned as well?\nWhen performing such high resolution analyses, and considering variables such as volume, these\ndifferent processes could have an impact. While I don’t think this would have too much of an\nimpact, it is important to at least document and be aware of it.\n--The molds measured by Pante were the second mold taken from each specimen.\nValitierra et al. 2023 used different methods of data collection and analysis than those that\nare applied here so their results are not comparable to this study. Their methods replicate\nthose of Maté-González et al., 2015, which used a low resolution structured light scanner\nand only data extracted from profiles. They also tested the impact of multiple moldings on\nmodern and not fossilized bones that have surfaces that are much more resistant to change\ndue to their lithified states. We agree with the reviewer that our molding process likely did\nnot impact our results, but we appreciate their research directed at determining the impact\nof molds on mark morphology.\nBeyond overlying trephic modifications, I also think it is very important to consider the impact\ndendritic etchings and root marks may have had on the samples, especially when considering the\nmethodology employed for quantitative BSM analysis and identification. For example, in\nVGr.0519, VGr.1557, VGr.1767, VGr.2275, VGr.2687 and VGr.2649 you can see overlapping\nbiochemical alterations that would surely affect the calculation of BSM volume and other\nparameters, wouldn’t these parameters now be unreliable for the identification of the traces? To\nwhat extent should the analysis therefore include these variables? What happens when you\nperform calculations on some of the samples including or excluding some of these parameters?\nDo the classifications change? It is very true that almost any methodology dedicated to the high-\nresolution analysis of BSMs would be conditioned by these alterations, however, calculations\nsuch as volume are more likely to be effected by such changes to the cortical surface.\n--The reviewer is correct that bioerosion can impact the measures collected here, which is\nwhy marks where bioerosive processes impacted depth were not included in the analysis.\nThe specimens cited above had multiple modifications present, and those that were altered\nwere excluded from our results.\nWith regards to data acquisition, much more information is needed regarding the microscopic\nanalysis, including the precise FOV and numerical aperture of the lens used and the windows\nscanned; variables that have been proven in the past by authors such as Calandra et al. (2019,\nDOI: 10.1038/s41598-019-42713-w) to effect the replicability of analyses.\n5\n--We have added the information for the objective used in our analysis. The 5x lens has a\nnumerical aperture of .15, a working distance of 23.5 mm, a field of view of 3400 um x 2837\num, a spatial sampling of 2.76 um, and an optical resolution of .93 um. However, we are not\nsure what the reviewer is referring to by “the window scanned”. It is important to note\nthat cut mark analysis is not equivalent to microwear analysis as the resolution required\nfor microwear is much higher and the measured variables are different.\nAlso, if I recall correctly, the S Neox has the option of obtaining 3D models using either Focus\nVariation or Confocal techniques, that would also condition the quality of the 3D model\nproduced, so the modality must be stated.\n--The S-Neox does not have a focus variation only mode, but it can combine focus variation\nwith confocal data or use confocal data alone. The results are not impacted by the different\ntechniques as the focus variation data is only incorporated to capture high slopes that are\nunmeasurable by the confocal method.\nWhy was only 5x magnification used? What colour light was used? My own experience with a\nsimilar microscope found 5x and 20x to produce different results for parameters such as\nroughness and volume when doing use-wear analysis. For this reason I’d also argue that\ncomparing data obtained from the S Neox and the Nanovea ST400 could definitely be considered\nproblematic for some of the quantitative variables included in the study (again, see the work of\nIvan Calandra). The authors state that the results should not be impacted because of “coarser\nscales of measurement” (which would probably answer my question of why was only 5x used),\nhowever this is yet to be tested, and when you calculate variables such as volume and roughness,\nthen multiple researchers, including myself, would in fact disagree with this statement.\n--We understand the reviewer’s concerns about the light and magnification, but it is\nimportant to keep in mind the scale of bone surface modification (BSM) analyses is much\nlarger than that for microwear. Microwear analyses only sample a small area of a tooth’s\nsurface, which allows the use of higher magnification objectives. BSM’s are many times\nlarger than this and high magnification measurements are not possible on this scale due to\nthe time required and limitations of the equipment used. For example, measuring a mark\nwith a 20x lens would produce a point cloud that would exceed the memory capabilities of\nthe S-Neox and would take hours instead of minutes to capture, if it were even possible.\nWe chose the 5x lens because it is the closest in resolution and xy sampling rate to the\npreviously used Nanovea st 400. Slight variations in measurements will occur in repeated\nmeasurements of marks even on the same device. This error range has been reported and\nexplored in Pante et al., 2017 (which is cited in the text currently) and is on the order of a\nfew microns for a variable like depth. Variations between measurements collected with\ndifferent color lights or between devices are on a nanometer scale and only relevant when\nconducting high-resolution surface texture analysis. For example, in our sample, the\nvariation between the mean values for depth between BSM are typically on the scale of tens\nto hundreds of microns so measurement variability between instruments and/or analysts\nequates to noise for our scales of analysis. This is even true for our measurements of Ra,\nwhich varies greatly between mark types. Ra measures the deviation of a surface from a\nmean line and the greatest deviation in a profile through a BSM is its overall depth.\nTherefore, our use of Ra is not in the traditional sense, which requires much higher\n6\nmagnification to measure. However, the variable as measured here remains useful in\ndiscriminating between mark types.\nInformation about how many non-measured points during the scans are also required\ninformation, so as to properly assess the quality of the scans used for the analysis.\n--Scans from molds generally result in 100% data capture. Occasionally missing data\npoints occur, but because they are isolated and outside of the mark they do not impact the\nmeasured variables. Every study that employs confocal microscopy fills in missing points\nand the quality of scans can be better judged by the provided image captures of the 3-D\nmodels, which are provided in the Supplemental figures.\nRegarding the post-processing of the 3D models obtained using the microscope, the author states\nthat they removed form from the data, yet no indication about what polynomial was used to do\nthis. I have to say that I think form removal for cut mark analysis is potentially very problematic,\nhowever I understand why it has been done considering some of the variables that have been\nmeasured. Form is removed in microwear studies because this directly has an impact on\nroughness parameters, and will also influence how easy it is to compute each traces volume,\nhowever when we are talking about the properties of the profile, form removal will have an\neffect on mark morphology.\n--These methods are described in Pante et al., 2017. Importantly, the mark is excluded\nfrom this transformation, which is detailed in Pante et al., 2017. Regardless, the\nquantitative method used effectively discriminates between the processed marks so\nremoving form is not problematic for our particular method and is necessary to capture\nthe variables of interest.\nWith regards to the statistical analysis of the data, and classification tasks performed, I also have\nsome concerns that need addressing. First of all, Quadratic Discriminant Analysis (QDA) should\nnot be used for ordination, as has been done for supplementary figure 22. While use of QDA\nalgorithms is perfectly acceptable for classification, it should not be used for ordination for a\nnumber of mathematical and epistemological reasons. QDA is a classification algorithm, so of\ncourse there will be separation between samples in the constructed and visualised feature spaces,\nbecause that’s what it is designed to do – force separation to exist. There is also a huge issue with\nthese types of constructed feature spaces due to issues with sample size and imbalance. While the\nindividual to variable ratio is acceptable for the dataset provided, the sample imbalance is large\nenough to be problematic, in particular between cut marks and percussion marks\n--While there is some literature that agrees with this it mostly relates to LDA and\ngeometric morphometrics and not QDA on variables related to size rather than shape. We\nemploy QDA here due to the unequal variance and size of the samples and we assess our\nmodel setting prior probabilities proportional to the occurrence of each mark to account\nfor the differences in sample size.\n(side note – why are these two compared in general, are there people who consider equifinality\nto exist between percussion and cut marks?).\n--Percussion marks inflicted by sharp edged anvils can resemble cut marks in our\nexperience. Short and shallow cut marks can also resemble microstriations associated with\n7\npercussion. The term equifinality can be misleading as it is a function of analyst experience\nwith control collections or the methods employed.\nWhat the authors should use instead, if they wish to visualise the multivariate differences\nbetween samples, is perform and visualise a PCA with a biplot, or simply present boxplots for\neach variable separately. If they want to use the QDA figure, then it must be used as a\ncomplementary representation of the data. Its also important to ensure that the feature space is\nconstructed using cross-validation, otherwise separation is much more likely to be spurious in\nnature. For an example, see https://doi.org/10.1007/s11692-020-09494-x - although this paper\nwas done for the purpose of discussing use of similar techniques in geometric morphometrics,\nthe algebraic properties of the ordination technique are very similar, while the input to the figures\nin this paper are not landmark data but instead standard multivariate distributions and are thus\ncompletely comparable to the present case study. Finally, an ordination technique like these ones\nshould not be presented without including in the axes labels the percentage of informative that is\ncaptured by each dimension, otherwise they are very hard to interpret.\n--We have added label percentages to the figure caption (Canonical 1 is 86.9, Canonical 2 is\n7.9%) and a PCA biplot as supplemental material.\nFor the evaluation of how well a classification algorithm works before being applied to\narchaeological data (presented in the first paragraph of page 22), the “accuracy” metric is often\nmisleading, especially in cases where the imbalance is notable. Instead the authors must present\nother metrics such as precision, recall, sensitivity or specificity, that are more informative of\nother parts of the confusion matrix\n--We have added an ROC curve as supplemental materials that reflects the requested\ninformation and shows the model is highly effective. Thanks for the suggestion! We will\nalso report accuracy as it is a standard that more scientists are familiar with.\nSome minor comments:\nIn the results section of the main article the authors state that a reduced number of bones\npresenting polishing is proof of limited transport. Polish does not necessarily imply movement,\nit’s the action of water-born abrasive particles altering the overall surface of the bone, this could\noccur with a low flowing course of water over a prolonged period of time without moving the\nbones at all. In the supplementary materials the authors mention rounding as well. I think it’s\nmore informative to simply present all of the water-based alterations to cortical surfaces in the\nmain text as well.\n--We have changed the term “transport” to “water-based alterations”\nIn the last paragraph of the taphonomy section in the results the authors mention that some\nspecimens have “more degraded surfaces”, it would be nice if the nature of this degradation be\nclarified – is it due to fluvial processes, for example?\n--We have added the clarification to specify that we mean root-etching and surface\nexfoliation.\nHave percussion marks or evidence of anthropogenic fracture patterns been identified in the site?\n8\n--Percussion marks were amongst the BSMs we were looking for, though no evidence of\nthem was found. We have updated the text to reflect this. We also recognize that breakage\npatterns are useful for taphonomic analysis, but due to our time constraints (we only had\n25 day total to examine 5,000+ specimens), we prioritize bone surface modifications in our\nanalysis. We acknowledge this and suggest it as a future study in our full taphonomic\nanalysis manuscript.\nOn line 227 the authors include some citations of papers that support the notion that sites exist\nwith identified cut marks but no associated lithic assemblages. Among these cited papers,\nthey include BK from Olduvai, however BK does contain a huge amount of lithics published as far\nback as the Leakeys, and is even mentioned in the abstract of the very paper they cite. I can\nassure you BK has lithics, while Mary Leakey even noted the presence of handaxes from this\nsite. You also cite Dikika, which is very problematic, and many would even argue that not only\nare there no lithics, but no cut marks either\n--The reviewer is corrected and we have removed the reference to BK from Olduvai and\nalso clarified our language to indicate that these are published reports of cut marks.\nIn summary, I once again wish to thank the authors for their work, I truly have no problem at all\nwith the validity of their finds, I just hope to ensure that the presentation of results are as\nrigorous, robust, and transparent, as possible. I have no problem disclosing my name as a\nreviewer - Lloyd A. Courtenay\n--We appreciate the positive comments and thoughtful suggestions. Thank you.\nReviewer #3 (Remarks to the Author):\nThe paper presented by Curran and colleagues, entitled \"Hominin presence in Eurasia by at least\n1.95 million years ago,\" presents clear evidence of ancient anthropogenic activity in the form of\nbutchery marks from several sites in the Olteț River Valley in Romania. The remains come\nmainly from a single Grăunceanu assemblage, although they also provide data from some other\npenecontemporaneous sites. A list of U-Pb dates establishes an age of 1.95 Ma. This work offers\nhighly relevant information for understanding human dispersal prior to 2 Ma and provides\npalaeoenvironmental data, through isotopic analysis, suggesting conditions conducive to the\ndevelopment of such human movements.\nThe authors present bone remains with unquestionable cut marks, which have been analyzed\nusing both quantitative and qualitative methods. Nevertheless, it is surprising that the authors\nstate that the fossils were examined under strong, low-angled light for evidence of bone surface\nmodifications. It would have been more appropriate to use microscopic observational methods\n(stereomicroscope, ESEM, digital microscopes, etc.). It seems plausible to suggest that this is the\nreason why the analysis method of Rodrigo et al., 2009, has been modified, with the aim of\navoiding the microscopic variables of the analysis. It is likely that had a microscopic approach\nbeen employed, a greater number of marks would have been identified. The marks attributed to\ncut marks are clearly visible to the naked eye, and therefore it is reasonable to hypothesize that a\nmicroscopic analysis would have identified a greater number of evidence, as well as helped to\nattribute some of the probable marks to high-confidence marks.\n9\n--We agree with the reviewer that the use of ESEM and/or stereomicroscopy would likely\nenable us to identify additional cut marks. Using macroscopic methods, hand-lenses, and\ndino-lites, we were able to identify 20 cut marks. Our team was constrained by time (and\nfunding) and used qualitative and quantitative methods to maximize our data collection on\nthe 5,000+ specimens in the time we had (~25 days total).\nThe authors propose that some of the striae are attributable to high-confidence cut marks, while\nothers may be associated with probable or unknown marks. Although they employ a quantitative\nmethodology, they utilize a similarity-by-association approach to ascertain the origin of the\nmarks as either anthropogenic or non-anthropogenic. It would be preferable for the authors to not\nuse association to attribute the marks, but rather to only apply the proposed analytical principles.\nIt is unclear why the teeth were excluded from the taphonomical analysis. If one of the objectives\nwas to identify butchering marks, the teeth may show evidence of cuts related to the cheek\ndefleshing and/or tongue removal. In addition, some of them could have been dated, obtaining\ndirect dates from remains modified by lithic tools.\n--We looked at any dentition that was still in the mandible or maxilla, but did not analyze\nisolated teeth, which we have clarified in the manuscript.\nIn the text (third paragraph of the Taphonomy section), the number of remains with high-\nconfidence cut marks is initially indicated as seven. However, later in the text, it is stated that\nthere are eight remains, and these are also the number of fossils that appear in Table S2.\nTherefore, it is necessary to correct the text.\n--Thank you for catching this- you are correct (there are seven high-confidence cut marks\nfrom VGr. and one from FM). We have made this clearer in the text.\nIn some instances, the observed marks appear to be clearly anthropic in origin. However, the\nanalysis suggests otherwise. This discrepancy may be a consequence of the macroscopic analysis\ndeveloped by the authors and the omission of certain microscopic variables that could have\nbiased the results. A good example is the metopodial VGr.0429. Upon examination of the\nimage (Fig. S12), it appears to be a clear anthropic disarticulation mark. However, the analysis suggests\nthat it is, in fact, a tooth mark. It is recommended that these probable cutmarks be reanalyzed\nunder microscopic tools (a stereomicroscope, at a minimum).\n--We agree with the reviewer that this mark is anthropic in origin. The mark was relatively\ndeep and on the edge of the condyle, which is the likely reason that the quantitative analysis\nassigned it to “tooth mark”. Our qualitative approach identified the mark as a cut mark\nand thus it is unlikely that further magnification would alter our conclusion. This\ncombination of analytical procedures (qualitative and quantitative) was used for all\nmarked specimens, thus we are confident in our assessment of the marks described here,\nrecognizing that it is conservative.\nVGr.1490 is described as a cutmark with a high probability and a high confidence, yet it appears\nboth in the text (page 6 of SI) and in table S2 as a probable cut mark. Could this be an error? In\nthe SI description, the authors even assign the activity carried out by hominins, tendon\nextraction.\n--VGr.1490 is identified as a “probable cut mark”. The quantitative analysis identified this\nas a cut mark (with high posterior probability), however given its curved trajectory, we\n10\nassigned the more conservative identification of “probable”, though we agree that is it most\nlikely anthropic in origin.\nWhy do the authors attribute the mark observed on specimen VGr.2052 (Fig. S18) to trampling if\nthe analysis indicates that it is a cut mark with a high probability? This suggests that the authors\nsometimes adhere to the result of the quantitative analysis, but at other times, the observation of\nthe researchers takes precedence.\n--We have edited our description of this mark and have clarified our final assessment\nprocedure, in which we did give slightly higher priority to the qualitative analysis done by\nhuman observers, who can implement a more holistic approach than the quantitative\nanalysis can alone.\nGiven the paucity of information on the stratigraphic context of these ancient collections, it is\nunderstandable that the authors refrain from offering interpretations regarding the type of\noccupation or the mode of access to the resources. However, it is crucial to clarify the authorship\nof the marks and their age.\n--We agree that these questions are important. However, our data cannot identify an author\nof the cut marks beyond hominin. We discussed the ages for the deposits in this manuscript.\nTables 1 and S7 are invaluable for the compilation of all ancient sites with evidence of hominids\nor anthropogenic activity (lithics and butchery marks) and the corresponding bibliography. It is\nrecommended to include additional bibliographical references concerning published extended\ntaphonomic studies of the cited sites:\n--We appreciate the reviewer’s compliment and have added these references.\n- For the Ain Boucherit site (Algeria):\nCáceres, I., Chelli Cheheb, R., van der Made, J., Harichane, Z., Boulaghraief, K., Sahnouni, M.,\n2023. Assessing the subsistence strategies of the earliest North African inhabitants: evidence\nfrom the Early Pleistocene site of Ain Boucherit (Algeria). Archaeological and Anthropological\nSciences 15, 87. https://doi.org/10.1007/s12520-023-01783-8\nFor the Sima del Elefante site (Spain):\nBlasco, R., Blain, H.A., Rosell, J., Díez, J.C., Huguet, R., Rodríguez, J., Arsuaga, J.L., Bermúdez\nDe Castro, J.M., Carbonell, E., 2011. Earliest evidence for human consumption of tortoises in the\nEuropean Early Pleistocene from Sima del Elefante, Sierra de Atapuerca, Spain. Journal of\nHuman Evolution 61, 503–509. https://doi.org/10.1016/j.jhevol.2011.06.002\nHuguet, R., Saladié, P., Cáceres, I., Díez, C., Rosell, J., Bennàsar, M., Blasco, R., Esteban-Nadal,\nM., Gabucio, M.J., Rodríguez-Hidalgo, A., Carbonell, E., 2013. Successful subsistence strategies\nof the first humans in south-western Europe. Quaternary International 295, 168–\n182. https://doi.org/10.1016/j.quaint.2012.11.015.\nThe evidence presented is of significant interest and should be published. The utilization of\nmicroscopy to analyze the attributes of the cut marks is recommended in order to identify\n11\nadditional marks, correctly assign probable ones, and enhance the images presented. It is also\nrecommended to add some high-resolution images of examples of trampling or carnivore\ntoothmarks.\n--We thank the reviewer for their insightful comments. We will be including further BSM\nimages in our other (detailed taphonomy) manuscript, which is currently under review at\nthe Journal of Human Evolution.\nReviewer #4 (Remarks to the Author):\nHominin presence in Eurasia by at least 1.95 million years ago, by Curran et al.\nThis is a very interesting work with some exciting results suggesting early evidence of human\npresence in East Europe by > 2 Ma, i.e. much older than thought until very recently. Please note\nthat the evaluation of the cut marks on bones and whether there are human made is beyond my\narea of expertise. I will mostly focus my review on the chronology and general context of the\nfossils.\nAs a general comment, I would say that the manuscript is overall well written and easy to follow.\nThe results and data presented apparently support the interpretation proposed by the authors.\nActually, the main issues are in my opinion not related to what is shown in the work, but rather\nto what is NOT shown. In that regard, the absence of a detailed contextualization of the fossils is\nmy major concern. I get that the fossils are coming from old excavations and that the various\nsites are most likely no longer accessible, but we have virtually no idea about their stratigraphic\nposition, and whether the fossils represent one single or various assemblages. In other words, the\nintegrity of the fossil assemblages is simply unknown. This major question mark about the\nintegrity of the fossil assemblage, and whether all fossils correspond to a single depositional\nevent, is actually not addressed in the present work, since the fossils dated with U-Pb are not\nthose that have cut marks. While one may assume that they have all the same age, we actually\ndon’t know, and this should be verified by running additional analyses on the cut-marked bones.\n--We have added further text to describe the depositional context, about which we have a\nfull manuscript under review at the Journal of Human Evolution. Our detailed taphonomic\nanalysis strongly suggests that these remains are from a single depositional context due to\nuniform surface appearance, lack of weathering, and overall consistent preservation. In\nsum, the remains are likely to have been deposited in a low-energy overbank flood. We\ntherefore have high confidence that these remains were deposited in a relatively short\nperiod of time and the U-Pb for Graunceanu also reflect the ages of the cut marked bones.\nOur U-Pb results for Fantana lui Mitilan and Valea Roscai also support the dating for\nGraunceanu since Fantana lui Mitilan has younger fauna and is dated to a younger age and\nValea Roscai has been compared to Graunceanu biochronologically and also has a similar\nage.\nThe same comment applies to the stable isotope analyses performed on the equid teeth. Direct U-\nseries dating should be performed in order to confirm its antiquity and make sure it can be\nclosely associated to the other fossils.\n12\n--Given the short time frame and lack of funding for further analyses, additional dating is\nnot possible at this time. Nevertheless, as discussed in the previous comment, the\ntaphonomic preservation of the VGr. assemblage is very uniform and thus, we consider it to\nbe a single deposit.\nIn that regard, I genuinely wonder whether Rare Earth Elements could also help demonstrate\nthat all fossils have the same origin and come from the same stratigraphic unit.\n--While we agree that rare earth elements (REE) may potentially be useful, in this context\nthe very large range of REE contents observed within a single tooth (as a result of variable\npostmortem uptake) we see no easy way that REE could be used in this way to distinguish\ndifferent stratigraphic sections.\nThe whole context is briefly summarized in Methods, but remains (too) vague. The authors\nmention that the location of various localities have been ‘reidentified’, therefore it should be\npossible to get further information about the context, and establish a proper stratigraphic column,\ndescription of the sediment, with some additional information about the formation of the fossil\nhorizon and the depositional environment for the 3 localities. A geological or geomorphological\nmap of the ORV with the stratigraphic correlations of the various sections would also be very\nuseful.\n--Though we have reidentified the site locations, the original excavations completely\nexhausted the fossil deposits. Nevertheless, the original reports for Graunceanu (the main\nfocus of this paper), describe the site as a single deposit (a “bone nest” or bed) and our\ntaphonomic analysis verifies this, as specimen surface condition was found to be extremely\nconsistent across the entire Graunceanu assemblage.\nThe lack of independent age control is also another major concern. There is some circular\nreasoning in using the biochronology (which is based on independently dated assemblages) to\nsupport the numerical ages obtained here. The uncertainty on biochronological inferences is\nusually minored and not really taken into account, and especially on large mammals. I doubt very\nmuch that the precision on the biochronological age is sufficient to undoubtedly position the site\naround 2.2-1.9 Ma, within a small 0.3 Ma age window. This would need further explanation with\na dedicated section in SM, and in particular about the robustness of the age control on the\nlocalities used as chronological tie points. Rodents may possibly help refining the chronology of\nthe assemblage (e.g. Maul & Markova QI, 2007), but seems like from Table S1 that there are too\nfew remains. As mentioned above, the absence of independent age control adds a non-negligible\nuncertainty to the dating results, especially when knowing that dental tissues behave as open\nsystem, and may show U leaching, which would artificially age the sample. I had a quick look at\nAndreescu et al. QI (2013), and I didn’t find the existing age control on the Tetoiu Fm very\nconvincing to securely confirm its chronological attribution to the Early Pleistocene.\n--Our research team has previously published a detailed biochronological assessment of the\nfauna at these fossil localities (Terhune et al., 2020); we have also published several\nadditional papers reviewing specific taxonomic groups (Werdelin et al., 2023; Croitor et al.,\n2024) that also discuss the biochronology of these fauna. Though many species from these\nlocalities are not useful from a chronological perspective, others (Sus strozzi, Bison\n(Eobison), Alces, Nyctereutes megamastoides, Vulpes alopecoides, Meles thorali, Pliocrocuta,\nCanis) serve as important biomarkers for the Early Pleistocene and consistently point\n13\ntoward Graunceanu being dated to the 2.2-1.9 Ma window we originally suggested in our\n2020 paper.\nThe figure below represents this biochronological assessment from Terhune et al. (2020)\n(though has been updated slightly to account for minor changes in taxonomy,\ncorresponding to Table S1). The gray box represents our approximate suggested time\nperiod of deposition for Graunceanu based on the fauna.\nThe U-Pb dates provide additional support to this estimate and are also consistent with\nbiochronology from the two other localities analyzed, Fantana lui Mitilan and Valea Roscai.\nWhile we very much wish we had micromammal remains with which to independently test\nthis dating, these fossils are lacking from the ORV localities. We count ourselves lucky that\nwe in fact have two relatively robust independent dating results for the ORV localities that\nare consistent with one another; as the reviewer likely already knows based on their last\ncomment for our Table S7 (see below), this is relatively rare. Indeed, many early Eurasian\nsites have been dated solely on the basis of biochronology or other relative dating\ntechniques like magnetostratigraphy (see updated Table S7, which now includes dating\nmethods for each site).\nother do show U-leaching (i.e., beyond equilibrium values). For example, this is the case for\nsample VGr.2256-2. My first reaction would be that this outcome may seriously question the\nreliability of the U-Pb age of 1.95 Ma obtained on that same tooth. I think this aspect would\ndeserve some additional discussion in the manuscript. From the information given in SM, I\nunderstand that both solution and laser ablation analyses were carried out on the teeth. Some\nsamples were microdrilled. I am not sure to understand why, and this could be clarified in the\ntext perhaps. It is also unclear to me how LA analyses were performed, either by drilling vertically within the\ndentine, or through mapping an area of the occlusal surface? If so, can we see the maps with the\nspatial distribution of the various ratios? Please add some extra explanations in SM. The exact\nlocation of the samples analyzed for U-series (given in Table S5) within the teeth could also be\ngiven. Fossil teeth are known to be highly heterogeneous in terms of spatial distribution of\nradioelement, and especially old teeth (e.g., Duval et al. QG, 2011) and it is unclear to me how\nthis uncertainty was evaluated. Some domains within the dentine or enamel may locally show\nvery distinct U-series data, which would bias the age results. It would be strongly advisable to\nsample other parts of each tooth to evaluate this variability and help identify more recent\nuranium uptake or leaching overprints.\n--We apologise for any confusion that may have arisen here. To clarify- all the U-Th\ndeterminations were performed on micro-drilled samples and all the U-Pb analyses were\nperformed via laser ablation. While we thought the former was clear from the existing\nsupplementary text, we have now added the words, ‘laser ablation’ to the main text U-Pb\nsection to further clarify.\nIt is also unclear to me how LA analyses were performed, either by drilling vertically within the\ndentine, or through mapping an area of the occlusal surface? If so, can we see the maps with the\nspatial distribution of the various ratios? Please add some extra explanations in SM. The exact\nlocation of the samples analyzed for U-series (given in Table S5) within the teeth could also be\ngiven. Fossil teeth are known to be highly heterogeneous in terms of spatial distribution of\nradioelement, and especially old teeth (e.g., Duval et al. QG, 2011) and it is unclear to me how\nthis uncertainty was evaluated. Some domains within the dentine or enamel may locally show\nvery distinct U-series data, which would bias the age results. It would be strongly advisable to\nsample other parts of each tooth to evaluate this variability and help identify more recent\nuranium uptake or leaching overprints.\n--Table S3 in the Supplementary information indicates that the laser sampling mode was\n‘static spot’ analysis. As a result there are no maps.\n--As we discussed in the SM, U mobility affecting U-Pb laser spot isochrons, and U mobility\naffecting 234U and 230Th ratios operate at different spatial and temporal scales in fossil\ndentine and are diagnosed in two different ways. We clarify this below but feel that the text\nwe currently provide in the supplementary discussion is sufficient.\n--The ultimate test for U-Pb isochron dating, particularly with such large spot populations\nis, “do we have a clean, clear U-Pb isochron?” If the answer is “yes,\" then we can infer\nwith confidence that an “early uptake” model is applicable, with relatively rapid post-\nburial U uptake from the surrounding environment and no subsequent U mobility on a\nscale that might significantly affect the equilibrium U-Pb age calculation. In this\nmanuscript we present U-Pb isochron data for a sub-set of the investigated teeth where this\nwas found to be the case.\n--We have used a simple reliability test for our 234U/238U measurements, by pairing them\nwith 230Th/238U measurements as for U-Th disequilibrium dating. Given the clearly\ngreater than 1 Ma age of these samples from their U-Pb isochrons, closed-system evolution\nfrom any plausible initial 234U/238U value must be so close to the infinite 238U-234U-\n230Th isochron so as to be analytically indistinguishable from it. Recent minor U gain or\n15\nloss (even where too minor to affect the U-Pb isochron age) will drive a sample away from\nthe infinite isochron, towards younger apparent U-Th ages (U gain) or to a clearly open-\nsystem position (U loss). Dentine 238U-234U-230Th measurements closer to the infinite\nisochron than their uncertainty indicate no recent U mobility at the point analysed and give\na reliable measurement of undisturbed 234U/238U for that tooth.\n--We measured 234U/238U on drilled dentine samples of ~ 2.5 mg and ~ 75 ug. The subset\nof the 75 ug analyses which had not gained or lost U gave a very similar mean 234U/238U\nto that of the ~2.5 mg samples, ~ 30x larger, indicating that recent dentine U mobility of\nthese teeth (which had all passed the U-Pb isochron test above) is fine-scale, and likely\ninternal. Any of the following would all lead to scatter in (and rejection of) the U-Pb\nisochron: greater U mobility, U mobility persisting over much of the burial time of the\ntooth, or significant heterogeneity of initial 234U/238U in the tooth. There is no expectation\nthat the tooth reported by Duval et al. (2011) would have passed this isochron test (had we\nbeen able to undertake this at the time the paper was written), as that tooth had clearly\nbeen subject to very significant recent U uptake unlike any of the teeth considered here and\nwould have been rejected from this study.\n--A successful (i.e., low scatter, low MSWD) U-Pb isochron indicates any U mobility is at a\nlevel which does not affect the U-Pb age. Minor, relatively recent U mobility can affect\nmeasured 234U/238U ratios so we have applied a simple test for this using the 238U-234U-\n230Th isochron space. Some of the 234U/238U sub-samples were found not to have been\nmeasurably affected by U mobility, and the tight clustering of their values allows us to\napply their mean (with sufficient uncertainty to cover their scatter) to all of our\ndisequilibrium corrected U-Pb isochrons.\nSpecific comments\n1. Main text, ‘The current best’: What does ‘best’ mean here? Please be more specific.\n--We have updated this wording.\n2. What is the definition of ‘ephemeral traces’ here? There several occurrences of this term\nthroughout the manuscript and I guess this should be explained somewhere.\n-- By ephemeral we are referring to the fleeting traces of hominins represented by either\nsmall or less-well documented localities that hint at a larger presence of hominins\nthroughout Eurasia in the Early Pleistocene, but that is not as robustly demonstrated as\nsites such as Dmanisi. We have tried to clarify this in the introduction.\n3. U-Pb dating, ‘2.01±0.20 Ma to 1.87±0.16’. Please indicate at least once in the main text\nwhether errors given are 1 or 2σ. ‘Ma’ missing after ‘1.87±0.16’\n--Values reported in the main text are ±2σ. This text has been updated and the missing Ma\nadded.\n4. U-Pb dating, ‘which is ubiquitous in fossil teeth and bone samples’: Do you mean at this\nspecific site? Please clarify the sentence. There are actually many works showing the timing of\n16\nthe U uptake may significantly vary within a given dental tissue or bone (e.g., Grün and Taylor\nAncient TL 1996; Grün, RM 2009; Duval et al., QG, 2011), between the various tissues of a\ngiven tooth, and between fossils from the same site/horizon.\n--Our statement refers to the fact that almost all fossil teeth and bone samples – from any\nsite – will experience some form of post mortem U-uptake. We are, of course, aware of the\nvariabilities in timing both between different sites and across individual samples but we do\nnot actually refer to timing at all here and so believe that our generic statement remains\naccurate.\n5. U-Pb dating, ‘Grăunceanu took place >1.95 Ma and potentially >2.01±0.20 Ma’: actually\nthese two ages mentioned here are virtually the same given the associated error. There is no need\nin my opinion to mention them both, one is sufficient.\n--We agree that from a statistical perspective these are the same, though we feel it’s\nimportant to emphasize that deposition of the site might have occurred well before 2.0 Ma.\nWe have edited this text slightly to reflect this nuance.\n6. U-Pb dating, ‘These results make these sites the most securely dated localities from Early\nPleistocene Europe’: the concept of ‘securely dated’ is very relative. The absence of independent\nage control and the presence of minimum ages on fossils do not make the sites better constrained\nthan others in my view.\n--We have updated this wording to now say “These U-Pb results, coupled with the\nbiochronology for these localities11-14, lend high confidence to our age of deposition for\nGrăunceanu of >1.95 Ma and potentially prior to 2.0 Ma, and indicate that Grăunceanu is\nthe earliest locality in Europe to show evidence for hominin activity.”\n7. Discussion, ‘Our paleoenvironmental data for Grăunceanu indicate a more open, arid\nenvironment with some nearby woodlands and water sources’ : do you mean more open and arid\nthan today? Please specify. Additionally, the estimated MAP for Grăunceanu is higher than\ncurrent conditions, which seems to be in contradiction with the statement above on a more arid\nenvironment.\n--We have changed our manuscript from reporting MAP to reporting seasonal data. That\nis, we estimated precipitation amounts for the driest (summer) and wettest (winter)\nmonths. Thus, the driest months would have been in summer, resulting in more arid\nconditions than today. We also removed “more” from that sentence.\n8. Discussion, ‘While it is clear that the hominin presence in Eurasia at this time was not\nextensive’: Why is it clear? And what do you mean by extensive? There is actually early\nevidence of human presence in Europe from Ukraine to Spain.\n--We have edited this text to say “While it is clear that the hominin presence in Eurasia at\nthis time was likely geographically and temporally discontinuous, the preponderance of\nephemeral traces for hominins in this region can no longer be ignored.”\n9. Fig. 3 caption: please specify whether errors are displayed at 1 or 2σ.\n--Values are reported as +/- Ma in Fig. 3.\n17\n10. Fig. S24: I am not sure to understand the interest of showing elemental concentrations in\ndental tissues? What kind of information does it bring?\n--We generated this image in our pilot studies and agree that it is not necessary for the\nmanuscript and have removed this image. Please note that this deletion does not show up as\na tracked change to allow for a more readable manuscript revision.\n11. Fig. S24: Metal concentration or elemental concentration?\n--Figure removed.\n12. Please would it be possible to show instead the maps of U concentrations, U234/U-238 and\napparent age for each of the dated samples?\n-- U-Th analyses were conducted by micro-drilling with very few sampling locations per\nsample and so maps cannot be constructed. The U-Pb dating method requires the use of all\nlaser ablation spots from a given sample in order to determine an age via the ‘isochron’\nconstruction. Once again, therefore, maps cannot be constructed using individual spots.\n13. Fog. S25: How was the present-day U234/U-238 measured?\n-- The determination of this parameter is explained in detail on page 8 of the\nSupplementary information.\n14. Table S1: Does the composition of the fossil assemblage come from previous studies (if so,\nplease cite the relevant references), or from the direct study of the specimens by tahe authors? If\nthe data given for this table are coming from a new taxonomic revision of the fossil assemblage,\nthen a dedicated section including a proper systematic paleontology study should be included.\n--This updated table incorporates our prior work describing the fauna from the Oltet River\nValley (Terhune et al., 2020, 2021; Croitor et al. 2024; Werdelin et al., 2023). These\nreferences are already provided in the main text where we cite this table, but we have also\nadded them to the caption for Table S1.\n15. Table S4: which U-234/U-238 values were used for the disequilibrium corrected age? This\nmay be clarified, as Table S5 shows that there are several values available for each sample.\n-- As noted above, page 8 of the Supplementary information contains a detailed explanation\nof how the value of 1.0039±0.0035, which is used for correction of all U-Pb ages, was\ndetermined.\n16. Table S5:\n- Please clarify whether MXF and UMF samples have all been measured through solution\nanalyses, or UMF samples have been analyzed using Laser Ablation?\n-- As indicated in the SM text and table, all 238U/234U measurements are by solution,\nper Hellstrom (2003). As noted in the table caption, MXF/UMF indicate the two sample\ndigestion procedures described in the SM text.\n- Where (=which lab) was Th-230/Th-232 ratio measured?\n-- The 230Th/232Th was missing a superscript “b” to indicate measurement at the\nUniversity of Melbourne after Hellstrom (2003).\n18\n- Not sure I understand well the U mobility column. Do ‘Closed’ mean ‘secular equilibrium’,\nand ‘U gain’ mean that the ratio are still below secular equilibrium?\n-- As noted in the table caption and SM text, and rephrased above, uranium mobility is\nindicated by measurable departure from the infinite age isochron. These samples cannot,\nhowever, be considered to be at secular equilibrium until 234U/238U is within its\nuncertainty of unity, which will occur hundreds of thousands of years into the future.\n- The apparent U-Th age can be added to this table. Out of curiosity, how ca we explain that the\nsame tooth may provide an apparent U-Th age of <500 ka and a U-Pb age of > 1.9 Ma. This\nwould probably deserve some additional comments in SM.\n-- We respectfully disagree with the reviewer here as we feel that the U-Th ages to this table\nwould be misleading. Where a U-Th age can be calculated, it simply indicates recent U gain\nfor the sub-sample analysed and has no meaning beyond that. Regarding the difference\nbetween calculated U-Pb age, and apparent U-Th age, we dedicated text in the SM to\ndiscussing this, also rephrased above. The U-Th age calculation, where it is possible to\nmake one, is not meaningful as it merely indicates recent micro-scale uranium gain.\n17. I have my doubts on the reliability of the age of many localities mentioned in that table. Can\nthe authors at least add one column mentioning the dating methods employed to constrain the\nage of the site (e.g., U-Pb, palaeomagnetism, biochronology, ESR, TCN, etc.)?\n--We have added the dating methods in a new column in Table S7.\nI have read the revised manuscript by Curran et al., titled “Hominin presence in Eurasia by at least 1.95\nmillion years ago.” As I mentioned in my previous review, the manuscript is generally well-written and of\nbroad interest. However, I am unable to comment on the robustness of the evidence for hominin\npresence or the U-Pb ages, as these techniques fall outside my expertise.\nMy main concern with the original version related to the interpretation of the stable isotope data. The\nauthors have provided a thorough and detailed response to my comments, and I am overall satisfied\nwith the revisions, which include:\na) Calculations of the oxygen isotope composition of meteoric precipitation 1.95 Ma using three\ndifferent transfer functions, along with an uncertainty estimate.\nb) A brief discussion justifying the conclusion of a higher winter precipitation component to aquifer\nrecharge relative to today.\nc) Estimates for precipitation amounts in the driest and wettest months, rather than an average for the\nwhole year.\nI would only recommend that the authors clarify in the text that the Kohn et al. equation was developed\nfor mean annual precipitation (MAP), and its validity for reconstructing seasonal precipitation is\nuntested. Additionally, while Kohn (2010) did not provide uncertainties for his equation, these are\ndiscussed in Kohn and McKay (2012, Paleo 3).\nOverall, I am satisfied with this revised version and recommend the manuscript for publication in Nature\nCommunications.\n--We thank the reviewer for their support of our paper. We have edited the methods and\nsupplemental text to include a discussion of the uncertainties in the seasonal range estimates, and\nadded the Kohn and McKay (2012) reference.\n_____________________\nResponse to Reviewer 2\nAuthor’s response to reviewer: The molds measured by Pante were the second mold taken from each\nspecimen.\nReply: Please include this in the paper, there is no reason why this very simple sentence cannot be\nincluded, thus responding to the reviewer’s (my) concerns.\n--We have added text to the methods section indicating that the second set of molds were employed\nin the quantitative analysis.\nanalysis than those that are applied here so their results are not comparable to this study. Their\nmethods replicate those of Maté-González et al., 2015, which used a low resolution structured light\nscanner and only data extracted from profiles.\nI debate: The authors don’t think there is a correlation between depth, width, and volume. This is a\npotentially mathematically flawed statement; please take a triangle, change its width and height, and\nsee how any minute change will also have an impact on the area. Regardless of whether the\nmeasurements taken are not precisely the same, if other authors note that there is a metric change to\ncut marks due to the moulding process, then one cannot completely disregard the study because slightly\ndifferent measurements are being taken by another study. Each of these metrics are inherently\ndependent on the global properties of a mark, even if measured on only a single profile, or the entire\nmark. Also, if other authors were able to observe changes in the metric properties of cut marks at such a\n‘low resolution’, then wouldn’t these changes be more pronounced at a higher resolution? Half the\ntime, lower resolution techniques are not able to differentiate between slight variations in variables, but\nhigh resolution techniques are, so if a low resolution technique is able to detect something, then it\nwould be reasonable to assume that the high resolution technique is more likely to detect such a\npattern. This can be simply acknowledged in the revised version of the text.\n--As noted above we have clarified in the text which set of molds were used in our analysis, and in\nresponse to the comments below we have clarified in the text that all bones are fully mineralized and\nunlikely to be affected by molding.\nAuthor’s response to reviewer: They also tested the impact of multiple moldings on modern and not\nfossilized bones that have surfaces that are much more resistant to change due to their lithified states.\nWe agree with the reviewer that our molding process likely did not impact our results, but we\nappreciate their research directed at determining the impact of molds on mark morphology.\nReviewer’s reply: I understand the reply, however the authors have not clarified anything in the text to\nsupport this. I have studied materials at 3 Ma that present no sign of mineralization at all, and are\nequally as fragile as a modern bone, while on the contrary I also currently study materials that are only\n20 ka that are completely mineralized. (1) Just because bones have fossilized, doesn’t make their cortical\nsurface more or less fragile, while (2) not all sites present the same degree of mineralization. The\nauthors haven’t mentioned at any point the mineralization of the bones in their paper. At this moment\nin time, we have to take the authors word that the moulding process has not damaged the bone. \nTransparency is needed and the authors should at least clarify some of these points. Also, if they are\ncertain that the two moulds are the same and no changes have occurred, then the first and the second\nmould of the same mark should produce the same metric results, could this not be tested, just in case?\n--We have added text in the methods section to clarify that all of the fossils analyzed as part of this\nstudy are heavily mineralized and that we did not observe any damage to the bones during the\n2\nmolding process. We therefore have high confidence that the fossils specimens were not damaged by\nthe molding process and that results from different molds are comparable.\nOverall reviewer’s reply – The authors reply with some valid points to me, the reviewer, but these points\nare still left unattended in the text, leaving all other readers of this body of work in the dark. They could\neasily clarify these points, in most cases with a single sentence. I repeat – if these concerns and\nquestions could have occurred to me, they can probably occur to somebody else. We need to be as\ntransparent as possible in our analyses, especially in the case of such extraordinary findings. Including\nthe sentence “The molds measured were the second mold taken from each specimen”, is a very simple\nstatement to include, with an additional acknowledgement about the potential risks, even if these risks\nmay not be entirely present. And if the authors are indeed certain that no changes have occurred, they\ncould easily check this by measuring the different moulds and comparing whether the measurement on\nmould 1, 2 and 3 are all the same, which would validate the point they are trying to make.\n--Please see comments to individual points above. At this point we feel that an additional analysis\nfocused on error related to molding individual marks is beyond the scope of this study. We also feel it\nlikely that any variation among molds is small in comparison to error associated with different\ntechniques or observers.\nAuthor’s response to reviewer: The reviewer is correct that bioerosion can impact the measures\ncollected here, which is why marks where bioerosive processes impacted depth were not included in the\nanalysis. The specimens cited above had multiple modifications present, and those that were altered\nwere excluded from our results.\nReviewer’s reply: have the authors clearly stated this in the text? It would be very easy to add a\nsentence that clarifies this. If this is stated in the text, please can they indicate where, because I cannot\nfind it.\n--We have added text to the methods section that states “in a few instances marks that were\nmolded correctly were subsequently removed from analysis due to bioerosion.”\nAuthor’s response to reviewer: We have added the information for the objective used in our analysis.\nReviewer’s reply – thank you for including this information\nAuthor’s response to reviewer: … For example, measuring a mark with a 20x lens would produce a point\ncloud that would exceed the memory capabilities of the S-Neox and would take hours instead of\nminutes to capture, if it were even possible…\nReviewer’s reply: I won’t go too much into this, but the idea that the point cloud exceeds the memory of\nthe computer is interesting... I, and other colleagues, study cut marks, tooth marks and percussion\nmarks at 20x, and have no problem at all with memory. Yes it takes longer, but science requires sacrifice,\nand while the work has been done on moulds, there is no issues with limited access to the samples\nwhich would imply that the authors could not have taken the time to scan all the marks at 20x. Maybe\nolder versions of the same microscope present storage issues that newer version of the S-Neox do not\npresent, if that is true then it might be necessary to state the year, version, or some indication about\nwhich edition of the microscope is being used, because I repeat, in our lab (University of Bordeaux), we\n3\nhave an S-Neox and have had no problem studying cut marks at 20x magnification.\n--The model of the S Neox scanner used in this study was manufactured in 2018 (now specified in the\ntext). Thus, the reviewer is correct that we used a slightly older version of the scanner that the\nreviewer also has access to; memory on this older version is limited in comparison to newer models.\nWe have added the manufacture date to the text. Given that our analyses with the 5x objective have\nclearly demonstrated anthropogenic alterations to the ORV remains, we feel that to reanalyze\neverything, including the sample of 898 experimentally made BSM, at a higher resolution would not\nprovide a return on the time investment that reanalysis would require. Our QDA model based on the\ndata collected with the 5x objective is shown through multiple statistical methods of assessment to be\nhighly effective at classifying BSM.\nAuthor’s response to reviewer: the quantitative method used effectively discriminates between the\nprocessed marks so removing form is not problematic for our particular method and is necessary to\ncapture the variables of interest.\nReviewer’s reply: I’ll have to take the author’s word for it, however they still have not stated anywhere\nwhich polynomial is used to remove the form from these scans. This is an incredibly simple statement to\ninclude in the methods – A polynomial of X was used to remove the form.\n--Methods for the removal of mark form are described in Pante et al., 2017 (which is already cited in\nthe methods text). Because the underlying shape of every bone is different, a different polynomial\nequation is used for each. We now specify that readers should see Pante et al., 2017 for further detail.\nAuthor’s response to reviewer: While there is some literature that agrees with [QDA should not be used\nfor ordination] it mostly relates to LDA and geometric morphometrics and not QDA on variables related\nto size rather than shape.\nReviewer’s reply: There is a lot to unpack, even in this short paragraph. First a point of debate:\n“variables related to size rather than shape” - Size and shape can be, and are often, inherently related -\nunless you provide evidence to say otherwise, this statement can be very problematic. Because the\nobjective of the paper is not to look at this question, I will brush past this statement, but I just would like\nto point out that care is needed here because this is not necessarily true, and in any other paper I would\nvery strongly disagree with such a statement. Regarding the much more important point, I think the\nauthor’s have missed the point with reference to the issues of using discriminant scores for ordination. I\nam not talking about geometric moephometrics (GMM), I am talking about mathematics. All the\nliterature I can find on QDA and the corresponding plots still relate them to the calculation of canonical\nvariates. In fact, QDA and LDA are built on some fundamentally identical algebraic calculations, that\nmeans they are completely related. Considering how QDA and LDA are heavily linked algebraically, this\nwould in turn imply that all issues related with LDA, which are related to CVA, are related in turn to QDA\nas well. I’m talking exclusively about algebra, not about GMM. The reason I mention GMM is simply\nbecause this is where the majority of research into these algebraic properties lie. However, I would like\nto point out that authors such as Cardini and Bookstein do not even use landmark based data in their\npapers about this topic, they use simple multivariate normal distributions, and thus still prove my point\n– QDA should not be used for ordination. Likewise, the paper by Albrecht (that I provide below) talks\nmostly about linear measurements, not landmarks… therefore it is directly relatable to the present\nstudy.\nIt is interesting that the scores in the plots are also called canonical something. From this perspective\nhow can you be so sure the plot is different from what has been criticised by other authors? Likewise, no\nindications or citations are included in the paper to clarify how calculations are performed, beyond the\nstatement that JMP was used. JMP is a private licensed software, which I do not have a licence for (in my\nopinion this in itself is an issue) so I cannot perform any tests to check for replicability, or even evaluate\nthe suitability and reliability of the algorithm in this specific context.\nIn sum; (1) if QDA is less sensitive to sample imabalance, i would like some references to support this, or\nsimulations to prove this. (2) You still need to clarify if cross validation is done, and if it isn’t performed,\nthen justification as to why must be provided, because in a lot of circumstances the exclusion of CV\nwould be considered completely problematic. (3) I am very happy to see a PCA is now included, which\nwould generally resolve a number of the issues I had with the QDA plot, however, there are a number of\nissues related to the PCA. (3a) The PCA provided is very very hard to read, and I cannot see where the\nfossil points fall. Have they even been included? Please improve the presentation of this figure so it is\nclear to read and see how the archaeological samples compare with the experimental ones. This is\nsomething you did with the QDA plot (you already know my opinion about how reliable this plot is),\nbut have refrained to do for the PCA; (3b) why have the confidence intervals been excluded for the PCA,\nbut included for the (problematic) QDA? Please standardise the presentation of your results. (3c) Having red\nand orange points on the PCA makes the samples really difficult to distinguish. Also, is it just me, or do\nthe colours not match up between the QDA and the PCA? Or even within the PCA at all. In the key there\nis blue, purple, green and purple, but in the PCA I cannot see blue, but instead see red, orange, green\nand purple. What is going on here?!\nIf the authors are still adamant that their use of QDA is appropriate, please provide justification. I\nprovide justification here to highlight why it is not appropriate (once again, before the authors claim\nthat literature originating from GMM is not relevant here, please read the methods sections of these\npapers to see that in most cases, the methods are not rooted in GMM analyses at all);\nAlbrecht, G. H. (1992). Assessing the affinities of fossils using canonical variates and generalized\ndistances. Journal of Human Evolution, 7(4), 49–69.\nBookstein, F. L. (2019). Pathologies of between-groups principal components analysis in geometric\nmorphometrics. Evolutionary Biology, 46, 271–302.\nCardini, A., & Polly, P. D. (2020). Cross-validated between Group PCA scatterplots: a solution to spurious\ngroup separation? Evolutionary Biology, 47, 85–95.\nCardini, A., O’Higgins, P., & Rohlf, F. J. (2019). Seeing distinct groups where there are none: spurious\npatterns from between-group PCA. Evolutionary Biology, 46, 303–316.\nMitteroecker, P., & Bookstein, F. (2011). Linear discrimination, ordination, and the visualization of\nselection gradients in modern morphometrics. Evolutionary Biology, 38, 100-114.\nRohlf, J. F. (2021). Why clusters and other patterns can seem to be found in analyses of high-\ndimensional data. Evolutionary Biology, 48, 1–16.\nIts incredibly important to point out that in each of these referenced papers, the authors may talk a lot\nabout GMM, however the data they are actually using to explore the algebraic and mathematical issues\npresented are simply multivariate distributions. This means the data used for these statistical and\ntheoretical experiments are completely analogous with the present author’s data. In particular I\nrecommend the paper by Cardini and Poly (2020) for proof on the importance of cross-validating their\ndiscriminant function plots. For the importance of cross-validation in classification I could recommend a\nlarge number of papers and books, however would simply state that the book by Hastie, Tibshirani and\nFriedman (2009) The Elements of Statistical Learning is more than a good starting point.\n--We appreciate the reviewer’s detailed response regarding QDA, LDA, CVA, and PCA. As we think is\ndemonstrated by the reviewer’s own text, there is considerable debate in the literature regarding the\nappropriate use of these methods, both in geometric morphometric analyses and beyond. Because of\nthis we feel it useful to provide both a CVA and PCA plot in our supplemental data, which will provide\nreaders with as much data as possible (as this reviewer advocates for throughout their response) with\nwhich to make informed interpretations of the data. Please also note that our CVA plot is an accurate\nrepresentation of the posterior probabilities that are used to assess the confidence of the QDA model\nin its classification of each mark. Marks that have posterior probabilities for multiple mark types\nappear in areas where the confidence intervals of these mark types overlap in the CVA plot. We do\nappreciate the reviewer’s comments about making sure that these two plots are consistent in their\npresentation; indeed we had thought the colors were comparable but some error must have occurred\nin the creation of the PCA plot that caused the problems noted. We have updated the figure and now\nhope that it meets the expectations of the reviewer.\nThe reviewer also questioned whether we used cross-validation, which was clearly reported in our\noriginal submission where we specified that a 25% validation set was extracted from the experimental\nsample of BSM to validate the model. We now call this a 25% cross-validation set for clarity.\nReply to the reviewer: We employ QDA here due to the unequal variance and size of the samples and\nwe assess our model setting prior probabilities proportional to the occurrence of each mark to account\nfor the differences in sample size.\nReviewer’s reply: Ok then please clarify this in the text, and also state what prior probabilities were\nused. This can be easily included in a single sentence.\n--Our existing text already states “Prior probabilities were set proportional to the occurrence of each\nmark type in the dataset to offset the disproportionate representations of each mark type in the\nexperimental sample.”\nReply to the reviewer: We have added an ROC curve as supplemental materials that reflects the\nrequested information and shows the model is highly effective. Thanks for the suggestion! We will also\nreport accuracy as it is a standard that more scientists are familiar with.\nReviewer’s reply: While the authors have already corrected the presentation of results to include more\nthan just accuracy (ROC curves are a great addition), I have to simply state, for the purpose of scientific\nclarification, that I strongly disagree with the term or concept that accuracy is a \"standard that more\nscientists are familiar with\" - I could not disagree more with the logic that just because everyone else\ndoes it means that it is right. For years it was general consensus that the world was flat, this doesn’t\nmean that these individuals are the future of science. For years nobody thought anything was smaller\nthan the atom, scientific research proved them wrong too (if anyone is paying attention they might\ncatch a Friends reference in here). Science evolves and so should we, just because \"scientists\"\n(this should actually be archaeologists/palaeoanthropologist/palaeontologists) are more familiar with\naccuracy, that doesn’t mean we shouldn't break the mould and develop our analyses to adopt correct\npractices. Moreover, if you actually consult the literature regarding classification in data science,\nthe authors would see that (mostly) non-archaeological scientists publish far more than simply accuracy.\nGiven this last point, I would argue that in data science reporting only accuracy is certainly not a\nstandard, which would make this statement incorrect on a number of different levels.\nJust to prove the point, and make sure that we go beyond the field of data science and into the realm of\narchaeology: imagine we have a new set of 100 marks; 80 cut, and 20 trampling. Imagine our algorithm\nclaims that all marks are cut marks. This would result in 80% accuracy, however, it would also be\nmisleading considering how the algorithm is actually allocating the label cut mark to everything (i.e.\nhasn't learnt to classify anything at all). Now, if we calculate sensitivity and specificity, at least one of\nthese values would be 0, successfully helping us detect that our 80% accuracy is misleading, and the\nalgorithm is in-fact awful, and should not be used in real case studies (EVER!). Now imagine that instead\nof archaeologists we are doctors and we have potentially killed 80 patients. Just because our patients\nhave been dead for 1.9 million years (wow, what a date!), doesn't mean we should not be just as\nrigorous and strict with our diagnostics. I apologise for rambling, but I find this point fundamental to\nclarify, especially if the authors find more great finds in the future and have different reviewers who are\nunfamiliar with this concept.\n--We have addressed these comments in our original revision by providing an ROC curve. We maintain\nthat classification accuracy is a valid measure for assessing the strength of discriminant models and\nalso a more familiar statistic to paleoanthropologists than an ROC curve.\n_____________________\nResponse to Reviewer 4\nThe authors have convincingly addressed most of my comments, and even when objecting, the answers\ngiven are reasonable and make sense. After going through the detailed rebuttal, it seems to me that the\nevidence provided on the three main critical aspects of work (i.e., human origin of the cut marks on\nbones; (ii) antiquity of the fossil material and (iii) palaeoenvironment reconstruction) are robust enough\nto support publication. I would especially like to thank the authors for dealing with the comments made\nby the reviewers in such a constructive way.\n-- We thank the reviewer for their continued support of this manuscript and appreciate their\nprofessional and constructive feedback.\nFirst, I must admit that I remain quite sceptical about the precision of the general chronostratigraphic\nframework of the ORV and Dacian basin. Yes, biochronology can give a rough chronological estimate of a\ngiven site, but this should be rather regarded as what it is, a rough estimation. Estimation because the\ninherent uncertainty typically associated to biochronology is usually minimized (or ignored). It may be\nrelated to (i) the taxonomic attribution of fossil specimens, which may vary depending on the specialist\ninvolved, (ii) the validity/robustness of large scale biostratigraphy correlations across a country or even a\ncontinent, which may often be questionable, and (iii) and the existing true independent age control\navailable at other sites used for comparison (i.e., based on methods like magnetostratigraphy and\nnumerical dating method), which is often limited and has an uncertainty on its own. This is why I\npersonally have difficulties to accept precise age constraints based on biochronology only (e.g., here:\n‘This formation extends from the base of the Pleistocene to as young as 1 Ma; Biochronological\nassessments indicate Grăunceanu is Late Villafranchian (~2.2-1.9 Ma) and is attributed to mammalian\nbiostratigraphic zones MN17/MmQ19-11). For example, in absence of numerical age results (and even\nthough, the existing age errors are usually not negligible) or magnetostratigraphic evidence indicating\nthe undisputed presence of Jaramillo, the minimum age constraint of 1 Ma proposed here is\nquestionable. Let’s say that an Early Pleistocene chronology may roughly be estimated for the Tetoiu Fm\nbased on biostratigraphic evidence. Following a similar reasoning, there are many examples of localities\nthat have been equally attributed to either the Middle, Late or Epivillrafranchian biochrons depending\non the authors. For example, I note that on the figure given in the rebuttal, the chronological window\nproposed for the site covers the latest Middle Villafranchian and earliest Late Villafranchian, whereas\nthe text mentions a much finer age window (‘Late Villafranchian (~2.2-1.9 Ma)’). Anyway, my point here\nis that I feel that the general chronostratigraphic context of the ORV and Dacian basin (and mainly the\nTetoiu Fm) would deserve a dedicated section in Supplementary Material, as biostratigraphy provides a\ncrucial support to the U-Pb dating results. This would help the reader unfamiliar with the area getting\nthe bigger picture. In particular, it is unclear how the various sites mentioned in the text do relate (e.g.,\nLa Pietriș, Valea Roșcăi, Fântâna lui Mitilan). It is mentioned that some of them are coeval, while other\nare younger/older, but no data are provided to support this. I understand that his information might be\nfound in various other papers, but it would be beneficial to the present work if a synthesis could be\ngiven in Supplementary Material. In particular the figure of the rebuttal showing the biochronological\nevaluation of the Grăunceanu should probably be added. It does provide useful information. The other\nlocalities form the ORV are that are mentioned in the text could be added to that figure too.\n--We agree completely that biochronological estimates of ages for fossil sites are challenging, and we\nhad the same concerns that the reviewer outlines here when we first described the biochronology of\nthese localities in our original publication (Terhune et al., 2020). This is in fact what spurred us on to\nseek the radiometric dates presented in this paper. Happily, our U-Pb dates and all prior\nbiochronological estimates, both by our team and previous researchers, are roughly consistent with\none another, and they are consistent with previous descriptions of the faunal horizons of the ORV\n(which include the sites discussed here) as described by Radulescu and Samson (1990). We agree with\nthe reviewer that more clarity of the positions and relationships of the sites and faunal horizons is\nvaluable and we have added a section in this regard to the Supplemental Materials. However, we\nhave decided not to include a detailed discussion of the biochronology of the fauna from the ORV\nsites, as it is our opinion that this would replicate previous work and take away from the focus of the\npaper, which is on the U-Pb dating and cut marks. We do agree with the reviewer that more precision\nregarding the upper limit of the Tetoiu Formation is necessary, and have also updated this in the main\ntext and supplemental to 1.3 Ma to be consistent with descriptions by Andreescu et al.\nFinally, I would like to add that, in my view, there is still an inherent general uncertainty on the context\nof the site and of the fossil assemblage. I do get the explanations of the authors on this matter, and I\naccept them, but the fact that the site is no longer accessible, and that most of the information related\nto the excavation is gone, naturally raises some questions on the integrity, extension and stratigraphy of\nthe fossil assemblage. This is complicated by the fact that various successive campaigns have been\nperformed over more than a decade, and we ignore how they relate. Moreover, dating analyses were\nnot performed on the same bones that have the cut marks, or on the equid teeth used for\npalaeoenvironment reconstruction, so there is also an inherent uncertainty on whether all these\nremains are truly associated, i.e., come from the same excavation area and the same stratigraphic unit.\nIn that regard, would it be possible to at least know from which excavation campaigns were taken all the\nspecimens analysed here for cut marks, dating and paleoenvironement analyses? That said, I get the\nauthors’ explanations and especially regarding the taphonomy suggesting a single depositional event,\nand no doubt they did their best to narrow down this uncertainty, but this is still indirect evidence that\ncannot fully compensate for the absence of a proper stratigraphy and spatial location of the remains.\n--We appreciate the reviewer’s comments here and understand where this confusion comes from. All\nof the materials examined for this study from Graunceanu come from the original excavations that\ntook place in the 1960s. Subsequent investigations have not revealed additional faunal accumulations\nat Graunceanu, since it seems likely that the entirety of the bone bed was exhausted in the original\nexcavations. Our taphonomic analyses (as well as original images and descriptions of these\nexcavations) strongly suggest that the faunal assemblage from Graunceanu stems from a single\ndepositional event. We very much wish we had a clearer picture of the internal stratigraphy of the\nsite, but in the absence of that we still have high confidence in the depositional integrity of the\nlocality. This is further supported by a continuous chain of custody of these remains from the original\nexcavators; one of those original excavators, Constantin Radulescu, was the PhD supervisor of Alex\nPetculescu, who is a coauthor on this study and currently manages the Bucharest collection. In the\nknown history of the collection there is no mention of anything other than the homogeneous site that\nwas excavated in the 1960s. Though the collection at present is split between two museums, it is very\nclear from markings on the bones that they were excavated at the same time and then subsequently\nsplit up in the intervening decades (no one seems to know why this split occurred).\n--We agree it would have been ideal to directly date the cut marked specimens; this certainly would\ngive us even higher confidence in the dates of any anthropogenic activity. However, methodologically\nit was not possible for this to happen, since the U-Pb dating requires dental samples and the vast\nmajority of the specimens exhibiting cut marks are from the postcrania. The one mandibular specimen\nwith cut marks (from Fantana lui Mitilan) is a nearly complete cervid mandible and dentition which is\nimportant both taxonomically and taphonomically, and therefore in our opinion should not be\ndestructively sampled. While we did destructively sample the horse dentition for isotopic study, this\nanalysis was conducted after the U-series dating campaign was concluded.\n--We have edited the methods in response to both of these comments and hope these edits\nsufficiently alleviate the reviewer’s concerns.\nSpecific comments\n1. L.74-75: The chronostratigraphy of the Tetoiu Fm described from Andreescu et al. (2013) is not as\nrobust and precise as suggested here. There is no independent age control available, and if I am not\nwrong the age constraints are only based on biostratigraphy. I’d rather see something like: ‘An Early\nPleistocene age may be considered for this Fm based on biostratigraphic evidence.’. I did not see\nevidence in Andreescu et al. (2013) for such a precise minimum age boundary of 1 Ma.\n--see comments above.\n2. L.104-106: Please indicate that this relative chronology is based on biostratigraphic evidence.\n--done\n3. Fig. 1: please check the spelling of ‘El-Kherba’.\n--Thanks for catching this error. The figure has been updated.\n4. L.196-198: I understand these ages only refer to Grăunceanu. This should probably be clarified.\n--Done\n5. L.199, ‘cluster tightly’: maybe a number (e.g., coefficient of variation, or a time range between the\nyoungest and oldest results) can be provided to support this statement.\n--A time range between the oldest and the youngest dates are already provided in the preceding\nsentence; we feel that adding anything else here might be more confusing than helpful.\n6. Fig. 3: please indicate in caption that age errors are 2σ.\n--Done\n10\n7. L. 231: VSMOW should probably be spelled out for its first occurrence in the text.\n--Done\n8. Table 1 and Supplementary Table S7: The various ages reported there are quite confusing, maybe this\nshould be clarified in caption. Sometimes an age range is given, with either 1 or 2 decimal digits. Not\nsure a 2 decimal digit precision is appropriate for most sites (Fântâna lui Mitilan even shows 3 digits). It\nis also unclear whether the age range is taken from the errors on the numerical ages, from the\ncombination of various ages, or simply correspond to various levels from the same site.\n--Thank you for this comment. It is very difficult to collate data from so many different localities and\nensure that the information presented is comparable. We have done our best in this regard but feel it\nis important to leave the ages of the sites as close to how they are reported in the original\npublications as possible. This means that some dates are reported with more precision than others,\nand others provide age ranges instead of specific numeric dates. We have tried to clarify this in the\ncaption for Supplemental Table 7. We have changed the date reported for Fantana lui Mitilan to more\nclosely reflect what is presented in our main text, however.\n- Not sure to understand how the age range of 2.3-2.0 Ma for Grăunceanu was obtained since the\nvarious U-Pb ages range from 1.87 to 2.01 Ma\n--Thank you for catching this error. We have edited this to now say >1.95 Ma.\n- Ain Hanech and El Kherba (Algeria) are usually regarded as coeval. Please see Duval et al. (2023) for the\nmost recent age estimation.\n--Thank you for pointing us to this publication. We have updated the ages for both sites in our tables.\n- Ain Boucherit (Algeria) is made of two levels, the lower level (Lw) dated to ~2.4 Ma and the upper level\n(Up) to ~1.9 Ma. Maybe these 2 different levels (usually referred as Ain Boucherit Lw and Ain Boucherit\nUp) should be differentiated, otherwise the age range given may suggest a significant uncertainty on the\nage of the site, which is not the case.\n--We have separated these two layers in both the table in the main text and Supplemental tables 7\nand 8.\n- Pirro Nord: please note the most recent dating study on that site (Duval et al., 2024). Although a\nbiochronological age of 1.3-1.6 Ma cannot be excluded, there is strong evidence suggesting that the site\nmay be much younger.\n--Thank you for pointing us to this publication. We have retained the original dates of 1.3-1.6 Ma, but\nhave also included a reference to this paper in the supplemental table.\n11\n- Sima del Elefante : maybe indicate the layers, e.g., Sima del Elefante (TE7-9). Methods also include BS.\n--We have opted not to include the level for this site since we generally do not provide level\ninformation for other sites (or at least not those with multiple distinct levels for which data are\nprovided). We have added BS as one of the dating methods.\n- Fuente Nueva-3: a rounded age of ‘1.2 Ma’ would probably be more appropriate given the associated\nuncertainty.\n--While we agree there is likely uncertainty here, we retain 1.19 Ma as the age since this is what is\nreported in the literature (see our comment above regarding the presentation of the age data for the\nfossil sites).\n- Barranco Leon: Methods also include MS and BS, in addition to ESR.\n--added\n9. L.324-330: for a matter of consistency, please use the same number of decimal digits for all\npercentages.\n--fixed\n10. L.326: Do the 3.32% include both Lower and Upper levels at Ain Boucherit? Please clarify.\n--That percentage did include both levels but we have now separated them out\n11. L.335: maybe you can refer to Supplementary Table S8 here.\n--added\n12. Ref 71 in SM: please note the name of the journal is ‘L'Anthropologie’.\n--fixed\n13. Supplementary Table S8: please indicate in caption the meaning of ‘Full Taph?’. Note that the age\nreported for Grăunceanu is not consistent with that in Table S7. Same comment applies for Dmanisi.\n--dates for Graunceanu and Dmanisi have been updated; thank you for catching those inconsistencies.\nFull Taph refers to whether a full taphonomic assessment of the entire assemblage for that site has\nbeen conducted and published. We have clarified that in the table.\nAgain it is a pleasure to be involved in the reviewing process of this paper, however I have rewritten\nthis review a number of times to make an effort to be as clear as possible. While some of the\ncomments made in response to my last review are debatable, they are not necessarily the focus of\nthe study, therefore I will leave these points for now, and this can hopefully lead to interesting\ndebates and lines of research in the future. Nevertheless, there is one point that has still not been\ncorrected to a sufficient extent, and is arguably one of the most important points I have made.\nWithout properly responding to this point, I do not feel comfortable accepting the metric analysis of\nthe taphonomic marks in their current state. Once again I have no problem in saying that the marks\nare cut marks, from a traditional perspective, but my issue is that when using quantitative\napproaches, a certain number of statistical and mathematical requirements have to be met in order\nfor the results to be considered robust, reliable, and not a misrepresentation of spurious\nseparations between samples.\nIn my first and second review, I stated very clearly that it is fundamental that QDA results are cross-\nvalidated. This is important for both the evaluation of the algorithms used for classification, and\ncan wildly change the quality of the ordination plots. When I brought it up again in the second round\nof reviews (and I provided a number of references to prove the importance of my point), the authors\nstate that it has already been included in their manuscript. JMP is a private licenced software, so I\nhave close to no way of testing some of the main issues, however, I have been able to obtain access\nto the non-pro version of JMP for a limited time in order to run some tests. Likewise, I have read in\ndetail the documentation of the software to try and find the information I requested. Finally, in this\nreview I will be providing R code to help clarify not only my points, but also provide a potential\nmeans of fixing some of the issues. I hope that the authors can see that the fact that I have taken\nthe time to provide an R script is not only for the purpose of proving the importance of cross-\nvalidation, but in the hope that the authors may also find this code useful. If the authors wish to use\ntheir own data, then I have provided comments within the code to clearly indicate where their data\nneeds to be inserted.\nThe authors state that in their original manuscript that they “clearly reported … that a 25%\nvalidation set was extracted from the experimental sample of BSM to validate the model. We now\ncall this a 25% cross-validation set for clarity”. I respectfully disagree, and I am not sure we are\ntalking about the same thing here. First, there are many types of cross-validation algorithms, which\none are you using (if you are using one at all)? In the case of classification, the most recommended\nform of performing cross-validation on a classification algorithm is k-fold cross-validation, with k =\n5, or more commonly k = 10. Did the authors do this? Next, there is a distinct difference between a\ntest set and a validation set, which is where I think some source of conclusion lies. Finally, how can\n25% of the data be used to cross-validate the plot? In general cross-validation of the discriminant\nfunction plots is slightly different.\nWhile a number of authors can be found throughout literature using the terms validation and test\nsets for different things, in general, the standard approach is to refer to the validation set as the\ndata used by the algorithm itself to define the discriminant function, and the test set is what we\nseparate and then use to evaluate how well this discriminant function is actually working (I think\nthis is what you refer to as validation, given that JMP uses this term in the Analysis/Predictive\nModelling/Make Validation Column menu). I suspect that the authors are calling the validation set\nwhat I (and many) would actually call the test set. This would imply that they are still not necessarily\nperforming cross-validation, which is why the simple statement “25% cross-validation” does not\nactually or fully clarify anything.\nThrough playing with JMP, and the additional code I provide in R, I am pretty certain that the authors\nare not performing the cross-validation required. JMP appears to be able to do everything I will\ndescribe, however in the non-pro version of the software I cannot actually test this fully.\nNevertheless, the documentation clearly states that K-fold cross validation is possible\n(https://www.jmp.com/support/help/en/18.0/index.shtml#page/jmp/kfold-and-holdback-\nvalidation.shtml#), which goes beyond simply separating a test set to evaluate model performance\nafterwards. This however requires a number of additional steps to be performed, however (see\nprevious link).\nFrom what I can tell by the table on this website, the algorithm is not doing k-fold cross validation\nautomatically, I can only assume that using a number of other functions they can get around this.\nOtherwise the word “No” under “K-fold cross-validation” and “Discriminant Models” does not bode\nwell for the use of this software for classification tasks using this algorithm. Likewise, I still do not\nsee clearly whether the cross-validation approach will be taken into account when computing the\ndiscriminant scores as well. This needs to be checked and clarified if the authors are adamant\nabout including a CV plot in their supplementary materials (which I have already suggested\navoiding all together).\nTo validate, and additionally clarify the points I am trying to make, I will describe the experiments I\nhave carried out in both R and JMP in continuation. If the authors are able to test that k-fold cross-\nvalidation is possible in JMP, and does indeed correct the plots, then this would be a good means of\nsolving all issues I am about to describe.\nIn the first example, I took the dataset from https://www.kaggle.com/datasets/muthuj7/weather-\ndataset?resource=download. This is a known benchmark dataset, and in R I can show that a cross-\nvalidated QDA is able to reach a 94% accuracy. I’d like to strongly point out however that (as I\nalready said in a previous review about how accuracy is a misleading metric), the QDA algorithm is\nunable to correctly classify any of the “null” class from this dataset, which leads to a kappa statistic\nof 0.65, i.e. the algorithm is not sufficient to correctly learn from this data. If I perform the same\nanalysis in JMP without performing k-fold cross-validation the results are completely different. The\naccuracy now drops to 45.9%, with a kappa statistic of 0.124. I think this is already evidence\nenough to show that cross-validation is not only a valuable means of boosting model performance,\nenhancing discriminant function robustness, but how misleading accuracy results can be if cross-\nvalidation is not performed. Confusion matrices of both these calculations are included as Figures\n1 and 2 in the file I provide alongside this review.\nTo provide sample sizes that are a little more analogous with archaeological datasets, however, I\nthen took a subset of 30 individuals per sample from this dataset. I used a seed in my R code to\nensure replicability of results, considering how random subsampling is performed. With this seed,\nyou will see that one of the observations included in the PCA (Fig. 3) is a potential outlier (it is in fact\nnot an outlier, as it does belong to the original distribution perfectly, a pretty cool example of how\nsampling strategies can make single extreme observations appear like outliers when they are in fact\nnot). In a cross-validated canonical variates plot (Fig. 4) the point is forced closer to the rest of the\nsample, however, a cross-validated plot (Fig. 5) shows that the point is still an outlier, which is\nmuch truer to the actual nature of the distribution. Either way, the two plots are different, showing\nthat cross-validation will produce a difference in plots and may force certain patterns to emerge\nthat may not be an accurate representation of the underlying data. To show how this effects cross-\nvalidation performance for classification, if we do not do cross-validation, JMP calculates an\naccuracy of 56.5%, with a Kappa of 0.347, while the cross-validated version of the algorithm\nproduces an accuracy of 56.5% and a Kappa of 0.3466. While the differences are minimal, they are\nthere.\nIn the second example I increase the complexity of the dataset by adding 100 new dimensions of\nnoise to the dataset, simply to provoke the ‘curse of dimensionality’ to occur, which again will lead\nto the spurious separation of samples mentioned by other authors. Once again note that without\ncross-validation, the plots are not comparable at all, and the cross-validated plot is much truer to\nthe original distribution as shown in the PCA. Likewise, the outlier from the first experiment is once\nagain still identified as such. In the non-cross-validated approach using JMP, the accuracy comes\nout at 43.5% but has a surprising kappa value of -0.049, indicating that the discriminant function is\nnot a good classifier at all. R, on the other hand, has the courtesy of refusing to perform the\nanalysis, as the sample size and complexity does not allow for the computation of a reliable\ndiscriminant function.\nIn the final example, I perform these experiments on purely synthetic data, that are designed to\nhave 0 differences between the three samples. Once again, spurious separation occurs on non-\ncross-validated results, but does not appear on the cross-validated results. This is proof of the\ntheorem that non cross-validated canonical variance plots can be incredibly difficult to interpret\nand are not always a necessarily reliable representation of the data. A non-cross-validated QDA in\nJMP reveals an accuracy of 47%, however once again R refuses to perform the computation.\nTo prove the fact that simply separating 75% of the data for the calculation of the discriminant\nfunction does not imply that cross-validation has been done, figures X, X, and X additionally show\nhow using the validation set and when excluding the validation set, spurious separation occurs\nregardless.\nIt is important to point out that in some of the above cases, only limited differences are seen\nbetween the results of k-fold cross-validated algorithms, and discriminant functions fit without\nperforming cross-validation. Regardless, the fact that even slight differences are there are\nimportant. In my previous reviews I provided a number of references to prove this. Seeing how this\nwas insufficient in stressing my point, here I have provided code, and statistical proof that this is\nstill an issue. Interestingly, R in a number of the above cases refused to calculate, or was unable to\ncalculate the discriminant function, while JMP succeeded. Nevertheless, I refer the authors to\nJMP’s documentation, where it clearly states that in some cases QDA can be problematic\n(https://www.jmp.com/support/help/en/18.0/index.shtml#page/jmp/discriminant-method.shtml#).\nThe authors need to ensure that they report robust and reliable quantitative results, and even if\nrepeating the statistics with a correct implementation of cross-validation produces the same\nresults, at least then we can be sure that your results are accurate.\nIn sum, for Sup. Fig. 22 please either; (1) cross-validate the calculation of the CV scores in the plot,\nor (2) completely delete this figure. In my opinion, option 2 is the easiest approach and does not put\ninto question the reliability of the study. For classification results, I insist that the authors must use\ncross-validation, preferably k-fold cross-validation with k = 10, as is standard procedure in most\nanalyses of this type. This latter point is the most important and needs to be rectified, or at least\nclarified.\nSigned – Lloyd A. Courtenay\nPS - I may be mistaken but I also recall an Excel or csv file being included in an earlier version of this\nmanuscript, which is now no longer available, if this was data it would have been nice to have\naccess to this data again to run the same tests on your actual data.\n--We appreciate the reviewer’s continued support of our manuscript and their thoughtful\ncomments on the statistical methods. It is clear from their comments that they feel strongly\nabout these methods, and we have done our best to address their concerns in as thorough a\nmanner as possible. We agree that some of the confusion here may have been related to\ndifferences in terminology, and we agree that a robust cross-validation of the model is critical\nfor ensuring proper interpretation of the results.\nTo that end, we have revised the text to make our analyses as clear as possible, and now\ninclude an additional author who performed updated analyses. As per the reviewer’s\nsuggestion we now present both leave-one-out and 10-fold cross-validation of the QDA\nresults. We have added a new Supplementary Table 3 that provides the confusion matrices,\nand have added the R code to run these analyses in the supplementary information (though\nnote this requires the raw experimental data, which are not yet published). Both cross-\nvalidation methods return very similar results (81.4% and 82%, respectively), demonstrating\nthe effectiveness of the QDA. In the course of these updates we noted a small error (n=3) in the\ncounts of bone surface modifications in the experimental dataset and have subsequently\nupdated Supplementary Figures 22 and 23, Supplementary Table 2, and all posterior\nprobabilities listed in the text. These changes are all very minor and do not change the\noutcomes of our interpretations; in fact, they strengthen them.\nAll changes made have either been tracked in the text or sections with major changes are\nhighlighted in yellow. Again, we appreciate the reviewer’s attention to our statistics and\nsupport of our work, and we hope these updates satisfy their concerns regarding the\nquantitative assessments of the bone surface modifications we present here.\nReviewer 2 attachments:\n# Example code for the cross-validation of discriminant plots and the\ntraining of classification algorithms\n# code by L.A.Courtenay, University of Bordeaux - [20 / 09 / 2024] to [24 /\n09 / 2024]\n# for queries or help contact ladc1995@gmail.com\n# Raw Example of the Weather History Dataset ------------\nlibrary(caret) # for cross-validation algorithms\nlibrary(MASS) # for QDA functions\nlibrary(Morpho) # the only package i know that offers cross-validated\ncanonical variate plots\nlibrary(ggplot2) # to visualise results\nset.seed(123456789)\n# !!!! Load here your own data, here i use a common known dataset just for\nan example !!!!!\n# the dataset i used can be downloaded from\nhttps://www.kaggle.com/datasets/muthuj7/weather-dataset?resource=download\ndataset <- read.table(\"weatherHistory.csv\", head = TRUE, sep = \",\")\ndataset <- dataset[,c(3, 4, 6, 7, 9, 11)] # remove some of the useless\nvariables that are not relevant for classification\n# !!!! Load here your own data, here i use a common known dataset just for\nan example !!!!!\n# extract numeric values and labels\nnumeric_values <- dataset[,2:6] # extract only numeric variables\nlabels <- as.factor(dataset[,1]) # extract the labels\n# perform a PCA\nnumeric_values_prima <- scale(numeric_values)\npca_results <- prcomp(numeric_values_prima)\nvariances <- pca_results$sdev ^ 2\nexplained_variance <- variances / sum(variances)\npc1_variance <- explained_variance[1] * 100\npc2_variance <- explained_variance[2] * 100\nplot_data <- data.frame(\npc1 = pca_results$x[,1],\npc2 = pca_results$x[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = pc1, y = pc2, colour = Sample)) +\ngeom_point(size = 1) +\nstat_ellipse(size = 1) +\nxlab(paste(\"PC1 (\", round(pc1_variance, 2), \"%)\")) +\nylab(paste(\"PC2 (\", round(pc2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n# because the dataset is huge, CVA will take an incredibly long time to\ncompute,\n# therefore for this example on the raw data I have refrained from doing so\nand\n# go straight onto calculating the benchmark accuracy using cross-validated\nQDA\n# seperating test data for evaluation\ntest_split <- 0.75\ntrain_x <- numeric_values[1:(ceiling(nrow(dataset) * test_split)),]\ntest_x <- numeric_values[ceiling(nrow(dataset) *\ntest_split):nrow(dataset),]\ntrain_y <- as.factor(labels[1:(ceiling(nrow(dataset) * test_split))])\ntest_y <- as.factor(labels[ceiling(nrow(dataset) *\ntest_split):nrow(dataset)])\ntrain_data <- data.frame(train_x, Sample = train_y)\ntest_data <- data.frame(test_x, Sample = test_y)\n# save the data so that it can be later imported into JMP\nwrite.table(\nrbind(data.frame(train_data, Validation = \"Train\"),\ndata.frame(test_data, Validation = \"Test\")\n),\n\"Full_JMP_data_for_QDA_experiments.csv\",\ncol.names = TRUE, row.names = FALSE, sep = \",\"\n)\n# cross-validated QDA model\ncross_validation_algorithm <- trainControl(method = \"cv\", number = 10)\nQDAmodel_CV <- train(\nSample ~ .,\ndata = train_data,\nmethod = \"qda\",\ntrControl = cross_validation_algorithm\n)\n# evaluate model\ntest_predictions <- predict(QDAmodel_CV, test_data)\nconfusionMatrix(test_predictions, test_y)\n#\n# Example 1 - Weather History dataset ----------------\nlibrary(caret) # for cross-validation algorithms\nlibrary(MASS) # for QDA functions\nlibrary(Morpho) # the only package i know that offers cross-validated\ncanonical variate plots\nlibrary(ggplot2) # to visualise results\nset.seed(123456789)\n# !!!! Load here your own data, here i use a common known dataset just for\nan example !!!!!\n# the dataset i used can be downloaded from\nhttps://www.kaggle.com/datasets/muthuj7/weather-dataset?resource=download\ndataset <- read.table(\"weatherHistory.csv\", head = TRUE, sep = \",\")\ndataset <- dataset[,c(3, 4, 6, 7, 9, 11)] # remove some of the useless\nvariables that are not relevant for classification\n# becase the dataset is huge i will also subsample 30 random examples from\neach of the classes\nsplit_data <- split(dataset, as.factor(dataset[,1]))\nsubsampled_data <- lapply(split_data, function(subset) {\nsubset[sample(nrow(subset), 30, replace = FALSE),]\n})\ndataset <- do.call(rbind, subsampled_data)\n# !!!! Load here your own data, here i use a common known dataset just for\nan example !!!!!\n# shuffle data\ndataset <- dataset[sample(nrow(dataset), replace = FALSE),]\n# extract numeric values and labels\nnumeric_values <- dataset[,2:6] # extract only numeric variables\nlabels <- as.factor(dataset[,1]) # extract the labels\n# perform a PCA\nnumeric_values_prima <- scale(numeric_values)\npca_results <- prcomp(numeric_values_prima)\nvariances <- pca_results$sdev ^ 2\nexplained_variance <- variances / sum(variances)\npc1_variance <- explained_variance[1] * 100\npc2_variance <- explained_variance[2] * 100\nplot_data <- data.frame(\npc1 = pca_results$x[,1],\npc2 = pca_results$x[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = pc1, y = pc2, colour = Sample)) +\ngeom_point(size = 4) +\nstat_ellipse(size = 1) +\nxlab(paste(\"PC1 (\", round(pc1_variance, 2), \"%)\")) +\nylab(paste(\"PC2 (\", round(pc2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n# non-cross-validated Canonical Variates plot\ncva_results <- CVA(numeric_values, labels, cv = FALSE)\ncv1_variance <- cva_results$Var[1,2][[1]]\ncv2_variance <- cva_results$Var[2,2][[1]]\nplot_data <- data.frame(\ncv1 = cva_results$CVscores[,1],\ncv2 = cva_results$CVscores[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = cv1, y = cv2, colour = Sample)) +\ngeom_point(size = 4) +\nstat_ellipse(size = 1) +\nxlab(paste(\"CV1 (\", round(cv1_variance, 2), \"%)\")) +\nylab(paste(\"CV2 (\", round(cv2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n# cross-validated Canoncial Variates plot\ncva_results_CV <- CVA(numeric_values, as.factor(labels), cv = TRUE)\ncv1_variance <- cva_results_CV$Var[1,2][[1]]\ncv2_variance <- cva_results_CV$Var[2,2][[1]]\nplot_data <- data.frame(\ncv1 = cva_results_CV$CVcv[,1],\ncv2 = cva_results_CV$CVcv[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = cv1, y = cv2, colour = Sample)) +\ngeom_point(size = 4) +\nstat_ellipse(size = 1) +\nxlab(paste(\"Cross Validated CV1 (\", round(cv1_variance, 2), \"%)\")) +\nylab(paste(\"Cross Validated CV2 (\", round(cv2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n# classification using QDA\n# seperating test data for evaluation\ntest_split <- 0.75\ntrain_x <- numeric_values[1:(ceiling(nrow(dataset) * test_split)),]\ntest_x <- numeric_values[ceiling(nrow(dataset) *\ntest_split):nrow(dataset),]\ntrain_y <- as.factor(labels[1:(ceiling(nrow(dataset) * test_split))])\ntest_y <- as.factor(labels[ceiling(nrow(dataset) *\ntest_split):nrow(dataset)])\ntrain_data <- data.frame(train_x, Sample = train_y)\ntest_data <- data.frame(test_x, Sample = test_y)\n# save the data so that it can be later imported into JMP\nwrite.table(\nrbind(data.frame(train_data, Validation = \"Train\"),\ndata.frame(test_data, Validation = \"Test\")\n),\n\"JMP_data_for_QDA_experiments.csv\",\ncol.names = TRUE, row.names = FALSE, sep = \",\"\n)\n# cross-validated QDA model\ncross_validation_algorithm <- trainControl(method = \"cv\", number = 10)\nQDAmodel_CV <- train(\nSample ~ .,\ndata = train_data,\nmethod = \"qda\",\ntrControl = cross_validation_algorithm\n)\n# evaluate model\ntest_predictions <- predict(QDAmodel_CV, test_data)\nconfusionMatrix(test_predictions, test_y)\n#\n# Example 2 - Weather History dataset ----------------\nlibrary(caret) # for cross-validation algorithms\nlibrary(MASS) # for QDA functions\nlibrary(Morpho) # the only package i know that offers cross-validated\ncanonical variate plots\nlibrary(ggplot2) # to visualise results\nset.seed(12345678)\n# !!!! Load here your own data, here i use a common known dataset just for\nan example !!!!!\n# the dataset i used can be downloaded from\nhttps://www.kaggle.com/datasets/muthuj7/weather-dataset?resource=download\ndataset <- read.table(\"weatherHistory.csv\", head = TRUE, sep = \",\")\ndataset <- dataset[,c(3, 4, 6, 7, 9, 11)] # remove some of the useless\nvariables that are not relevant for classification\n# becase the dataset is huge i will also subsample 30 random examples from\neach of the classes\nsplit_data <- split(dataset, as.factor(dataset[,1]))\nsubsampled_data <- lapply(split_data, function(subset) {\nsubset[sample(nrow(subset), 30, replace = FALSE),]\n})\ndataset <- do.call(rbind, subsampled_data)\n# add some artificial noise to make the dataset more interesting at higher\ndimensions\nfor (i in 1:100) {\ndataset <- cbind(dataset, rnorm(nrow(dataset), runif(1, 0, 10), runif(1,\n0, 10)))\n}\n# !!!! Load here your own data, here i use a common known dataset just for\nan example !!!!!\n# shuffle data\ndataset <- dataset[sample(nrow(dataset), replace = FALSE),]\n# extract numeric values and labels\nnumeric_values <- dataset[,2:ncol(dataset)] # extract only numeric\nvariables\nlabels <- as.factor(dataset[,1]) # extract the labels\n# perform a PCA\nnumeric_values_prima <- scale(numeric_values)\npca_results <- prcomp(numeric_values_prima)\nvariances <- pca_results$sdev ^ 2\nexplained_variance <- variances / sum(variances)\npc1_variance <- explained_variance[1] * 100\npc2_variance <- explained_variance[2] * 100\nplot_data <- data.frame(\npc1 = pca_results$x[,1],\npc2 = pca_results$x[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = pc1, y = pc2, colour = Sample)) +\ngeom_point(size = 4) +\nstat_ellipse(size = 1) +\nxlab(paste(\"PC1 (\", round(pc1_variance, 2), \"%)\")) +\nylab(paste(\"PC2 (\", round(pc2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n# non-cross-validated Canonical Variates plot\ncva_results <- CVA(numeric_values, labels, cv = FALSE)\ncv1_variance <- cva_results$Var[1,2][[1]]\ncv2_variance <- cva_results$Var[2,2][[1]]\nplot_data <- data.frame(\ncv1 = cva_results$CVscores[,1],\ncv2 = cva_results$CVscores[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = cv1, y = cv2, colour = Sample)) +\ngeom_point(size = 4) +\nstat_ellipse(size = 1) +\nxlab(paste(\"CV1 (\", round(cv1_variance, 2), \"%)\")) +\nylab(paste(\"CV2 (\", round(cv2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n# cross-validated Canoncial Variates plot\ncva_results_CV <- CVA(numeric_values, as.factor(labels), cv = TRUE)\ncv1_variance <- cva_results_CV$Var[1,2][[1]]\ncv2_variance <- cva_results_CV$Var[2,2][[1]]\nplot_data <- data.frame(\ncv1 = cva_results_CV$CVcv[,1],\ncv2 = cva_results_CV$CVcv[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = cv1, y = cv2, colour = Sample)) +\ngeom_point(size = 4) +\nstat_ellipse(size = 1) +\nxlab(paste(\"Cross Validated CV1 (\", round(cv1_variance, 2), \"%)\")) +\nylab(paste(\"Cross Validated CV2 (\", round(cv2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n#\n# seperating test data for evaluation\ntest_split <- 0.75\ntrain_x <- numeric_values[1:(ceiling(nrow(dataset) * test_split)),]\ntest_x <- numeric_values[ceiling(nrow(dataset) *\ntest_split):nrow(dataset),]\ntrain_y <- as.factor(labels[1:(ceiling(nrow(dataset) * test_split))])\ntest_y <- as.factor(labels[ceiling(nrow(dataset) *\ntest_split):nrow(dataset)])\ntrain_data <- data.frame(train_x, Sample = train_y)\ntest_data <- data.frame(test_x, Sample = test_y)\n# save the data so that it can be later imported into JMP\nwrite.table(\nrbind(data.frame(train_data, Validation = \"Train\"),\ndata.frame(test_data, Validation = \"Test\")\n),\n\"JMP_noisy_data_for_QDA_experiments.csv\",\ncol.names = TRUE, row.names = FALSE, sep = \",\"\n)\n# cross-validated QDA model\ncross_validation_algorithm <- trainControl(method = \"repeatedcv\", repeats =\n10)\nQDAmodel_CV <- train(\nSample ~ .,\ndata = train_data,\nmethod = \"qda\",\ntrControl = cross_validation_algorithm\n)\n# note that the algorithm gives an error because it is unable to compute a\ndiscriminant function\n#\n# Synthetic Data ----------------\nlibrary(caret) # for cross-validation algorithms\nlibrary(MASS) # for QDA functions\nlibrary(Morpho) # the only package i know that offers cross-validated\ncanonical variate plots\nlibrary(ggplot2) # to visualise results\nlibrary(abind) # for tensor manipulation operations\nset.seed(123456789)\n# !!!! Load here your own data, here i use a common known dataset just for\nan example !!!!!\n# for this example lets use purely synthetic data\ndataset <- array(numeric(), dim = c(90, 0))\nfor (i in 1:100) {\ndataset <- abind(dataset, rnorm(nrow(dataset), runif(1, 0, 10), runif(1,\n0, 10)),\nalong = 2)\n}\ndataset <- data.frame(labels = c(rep(1, 30), rep(2, 30), rep(3, 30)),\ndataset)\n# !!!! Load here your own data, here i use a common known dataset just for\nan example !!!!!\n# shuffle data\ndataset <- dataset[sample(nrow(dataset), replace = FALSE),]\n# extract numeric values and labels\nnumeric_values <- dataset[,2:ncol(dataset)] # extract only numeric\nvariables\nlabels <- as.factor(dataset[,1]) # extract the labels\n# perform a PCA\nnumeric_values_prima <- scale(numeric_values)\npca_results <- prcomp(numeric_values_prima)\nvariances <- pca_results$sdev ^ 2\nexplained_variance <- variances / sum(variances)\npc1_variance <- explained_variance[1] * 100\npc2_variance <- explained_variance[2] * 100\nplot_data <- data.frame(\npc1 = pca_results$x[,1],\npc2 = pca_results$x[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = pc1, y = pc2, colour = Sample)) +\ngeom_point(size = 4) +\nstat_ellipse(size = 1) +\nxlab(paste(\"PC1 (\", round(pc1_variance, 2), \"%)\")) +\nylab(paste(\"PC2 (\", round(pc2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n# non-cross-validated Canonical Variates plot\ncva_results <- CVA(numeric_values, labels, cv = FALSE)\ncv1_variance <- cva_results$Var[1,2][[1]]\ncv2_variance <- cva_results$Var[2,2][[1]]\nplot_data <- data.frame(\ncv1 = cva_results$CVscores[,1],\ncv2 = cva_results$CVscores[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = cv1, y = cv2, colour = Sample)) +\ngeom_point(size = 4) +\nstat_ellipse(size = 1) +\nxlab(paste(\"CV1 (\", round(cv1_variance, 2), \"%)\")) +\nylab(paste(\"CV2 (\", round(cv2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n# cross-validated Canoncial Variates plot\ncva_results_CV <- CVA(numeric_values, as.factor(labels), cv = TRUE)\ncv1_variance <- cva_results_CV$Var[1,2][[1]]\ncv2_variance <- cva_results_CV$Var[2,2][[1]]\nplot_data <- data.frame(\ncv1 = cva_results_CV$CVcv[,1],\ncv2 = cva_results_CV$CVcv[,2],\nSample = labels\n)\nggplot(data = plot_data,\naes(x = cv1, y = cv2, colour = Sample)) +\ngeom_point(size = 4) +\nstat_ellipse(size = 1) +\nxlab(paste(\"Cross Validated CV1 (\", round(cv1_variance, 2), \"%)\")) +\nylab(paste(\"Cross Validated CV2 (\", round(cv2_variance, 2), \"%)\")) +\ntheme_bw() +\ncoord_fixed()\n#\n# seperating test data for evaluation\ntest_split <- 0.75\ntrain_x <- numeric_values[1:(ceiling(nrow(dataset) * test_split)),]\ntest_x <- numeric_values[ceiling(nrow(dataset) *\ntest_split):nrow(dataset),]\ntrain_y <- as.factor(labels[1:(ceiling(nrow(dataset) * test_split))])\ntest_y <- as.factor(labels[ceiling(nrow(dataset) *\ntest_split):nrow(dataset)])\ntrain_data <- data.frame(train_x, Sample = train_y)\ntest_data <- data.frame(test_x, Sample = test_y)\n# export for JMP\nwrite.table(\nrbind(data.frame(train_data, Validation = \"Train\"),\ndata.frame(test_data, Validation = \"Test\")\n),\n\"JMP_noisy_data_for_QDA_experiments.csv\",\ncol.names = TRUE, row.names = FALSE, sep = \",\"\n)\n# note that R will not calculate the QDA for this data, therefore has been\nexcluded here\n# Classifying archaeological samples -------------------------\n# if the user wishes to calculate the labels for new specimens, such as the\narchaeological traces\n# you need to load the data using the read.table function, or have them\nseparated at the begining\n# then use the predict() function, e.g. predict(QDAmodel,\narchaeological_data),\n# to predict the new labels of each of the observations.\n# Prediction probabilities can be obtained by including the parameter type\n= \"prob\"\n# in the predict() function, e.g. predict(QDAmodel,\narchaeological_data, type = \"prob\")\n# for ROC curves, you can use the pROC or ROCR packages, however ROC curves\nare only\n# useful for comparing binary classification tasks. For categorical\nclassification,\n# you would need to calculate a curve for each of the pairs of labels. the\nroc()\n# function can be used to calculate the curve, and the auc() funciton will\nthen calculate\n# the area under the ROC curve. If sample imbalance is present, however,\nplease note\n# that ROC may not be the best approach for model evaluation. use help(roc)\nto\n# consult R's official documentation for this function.\nFigure 1 – A PCA of the dataset used for some of the experiments described in this review. Note that\nin a previous review the authors stated that because they do not use GMM data the comments made\nabout the reliability of cross-validation were not applicable to their own case study. In this\nexperiment I have not used, in any way, shape, or form, morphological data. The variables included\nare variables relating to the weather, humidity, and speed of the wind.\nFigure 1 – A screen shot of the confusion matrix produced by JMP when used to classify the dataset\nused for the first couple of experiments. For this confusion matrix, QDA was trained on 75% of the\ndata, and the confusion matrix was computed using the remaining 25% of the data for testing.\nFig. 2 – A screen shot of the confusion matrix produced by R when using a k-fold cross-validated (k =\n10) QDA algorithm. Once again the confusion matrix was computed using 25% of the data for testing.\nFig. 3 – A PCA presenting a subset of the original weather dataset, used for the first primary\nexperiment. If you consult Figure 1, you will see that the apparent outlier is not, in fact, an outlier,\nyet simply a product of sampling bias.\nFig. 4 – A Canonical Variance plot describing the separation of samples without the use of cross-\nvalidation.\nFig. 5 – A cross-validated Canonical Variance plot describing the separation of samples – note that\nthe outlier is still represented as such, implying that the plot is a more realistic representation of the\nactual underlying structure and nature of the data\nFig. 6 – A simple Canonical Variance plot describing the separation of the same dataset from figures\n3, 4, and 5, however 100 dimensions of noise have been added to provoke the spurious separation of\nsamples, proving that non-cross-validated plots are not necessarily a robust or reliable\nrepresentation of the data.\nFig. 7 - A cross-validated Canonical Variance plot of the same data from Figure 6. Note that the\ndistribution of samples is completely different and not at all comparable, even though they were\ncomputed on the exact same dataset.\nFig. 8 - A simple Canonical Variance plot describing the separation of a purely synthetic 100\ndimensional dataset where absolutely no differences between samples is supposed to exist. All three\nsamples come from the exact same probability distribution. This is theoretical proof that non-cross-\nvalidated plots are not a reliable representation of these distributions, as they are supposed to all be\nthe same. While this is product of an inflated variable to individuo ratio, it is an exaggerated means\nof proving the Marchenko-Pastur theorem.\nFig. 9 - A cross-validated canonical variance plot of the data presented in figure 7. Note that the plots\nare not at all comparable.\nFig. 10 – An example of calculating CV plots on the data from Figure 7 without performing the\nvalidation split. This proves that JMP is not necessarily cross-validating the computation of canonical\nvariance scores, and thus the distributions presented in the plot are not reliable.\nFig. 10 – An example of calculating CV plots on the data from Figure 7 while using the 75% training\nsplit. The only reason tis plot and the plot from Figure 10 are different is simply because the number\nof data points is lower, considering how 25% of the data has been excluded for “validation” (or\ntesting). This is further proof that simply using the nalysis/Predictive Modelling/Make Validation\nColumn does not solve any of the issues or points I have been trying to make since my very first\nreview letter."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have read the revised manuscript by Curran et al., titled “Hominin presence in Eurasia by at least 1.95 million years ago.”\nAs I mentioned in my previous review, the manuscript is generally well-written and of broad interest. However, I am unable to\ncomment on the robustness of the evidence for hominin presence or the U-Pb ages, as these techniques fall outside my\nexpertise.\nMy main concern with the original version related to the interpretation of the stable isotope data. The authors have provided\na thorough and detailed response to my comments, and I am overall satisfied with the revisions, which include:\na) Calculations of the oxygen isotope composition of meteoric precipitation 1.95 Ma using three different transfer functions,\nalong with an uncertainty estimate.\nb) A brief discussion justifying the conclusion of a higher winter precipitation component to aquifer recharge relative to today.\nc) Estimates for precipitation amounts in the driest and wettest months, rather than an average for the whole year.\nI would only recommend that the authors clarify in the text that the Kohn et al. equation was developed for mean annual precipitation (MAP), and its validity for reconstructing seasonal precipitation is untested. Additionally, while Kohn (2010) did not provide uncertainties for his equation, these are discussed in Kohn and McKay (2012, Paleo 3).\nOverall, I am satisfied with this revised version and recommend the manuscript for publication in Nature Communications.", "reviewer_2": "I am grateful to once again have the chance to review the paper by Curran and colleagues, and am happy to see that they\nhave replied in a constructive way to each of my comments and concerns, clarifying or rectifying the majority of the points I\nraised. There are some points, however, that have not fully been touched upon, while others that the authors have well\njustified in their letter to me, the reviewer, but should still be clarified in the text. I would just like to point out that if these\nquestions have occurred to me, they can easily occur to somebody else, and our objective should therefore be to have a text\nthat is as transparent and informative as possible, so as to avoid any further scepticism or questions. I reply here to some\npoints that I still think are worth debating, but are not entirely an issue if they are not confronted in the text, while I insist on\nsome other points that I feel are fundamental and must still be clarified.\n\nAuthor’s response to reviewer: The molds measured by Pante were the second mold taken from each specimen.\nReply: Please include this in the paper, there is no reason why this very simple sentence cannot be included, thus\nresponding to the reviewer’s (my) concerns.\n\nAuthor’s response to reviewer: Valitierra et al. 2023 used different methods of data collection and analysis than those that\nare applied here so their results are not comparable to this study. Their methods replicate those of Maté-González et al.,\n2015, which used a low resolution structured light scanner and only data extracted from profiles.\n\nI debate: The authors don’t think there is a correlation between depth, width, and volume. This is a potentially\nmathematically flawed statement; please take a triangle, change its width and height, and see how any minute change will\nalso have an impact on the area. Regardless of whether the measurements taken are not precisely the same, if other authors\nnote that there is a metric change to cut marks due to the moulding process, then one cannot completely disregard the study\nbecause slightly different measurements are being taken by another study. Each of these metrics are inherently\ndependent on the global properties of a mark, even if measured on only a single profile, or the entire mark. Also, if other authors were\nable to observe changes in the metric properties of cut marks at such a ‘low resolution’, then wouldn’t these changes be\nmore pronounced at a higher resolution? Half the time, lower resolution techniques are not able to differentiate between\nslight variations in variables, but high resolution techniques are, so if a low resolution technique is able to detect something,\nthen it would be reasonable to assume that the high resolution technique is more likely to detect such a pattern. This can be\nsimply acknowledged in the revised version of the text.\n\nAuthor’s response to reviewer: They also tested the impact of multiple moldings on modern and not fossilized bones that\nhave surfaces that are much more resistant to change due to their lithified states. We agree with the reviewer that our\nmolding process likely did not impact our results, but we appreciate their research directed at determining the impact of\nmolds on mark morphology.\nReviewer’s reply: I understand the reply, however the authors have not clarified anything in the text to support this. I have\nstudied materials at 3 Ma that present no sign of mineralization at all, and are equally as fragile as a modern bone, while on\nthe contrary I also currently study materials that are only 20 ka that are completely mineralized. (1) Just because bones have\nfossilized, doesn’t make their cortical surface more or less fragile, while (2) not all sites present the same degree of\nmineralization. The authors haven’t mentioned at any point the mineralization of the bones in their paper. At this moment\nin time, we have to take the authors word that the moulding process has not damaged the bone. Transparency is needed and\nthe authors should at least clarify some of these points. Also, if they are certain that the two moulds are the same and no\nchanges have occurred, then the first and the second mould of the same mark should produce the same metric results, could\nthis not be tested, just in case?\n\nOverall reviewer’s reply – The authors reply with some valid points to me, the reviewer, but these points are still left\nunattended in the text, leaving all other readers of this body of work in the dark. They could easily clarify these points, in\nmost cases with a single sentence. I repeat – if these concerns and questions could have occurred to me, they can probably\noccur to somebody else. We need to be as transparent as possible in our analyses, especially in the case of such\nextraordinary findings. Including the sentence “The molds measured were the second mold taken from each specimen”,\nis a very simple statement to include, with an additional acknowledgement about the potential risks, even if these risks may not\nbe entirely present. And if the authors are indeed certain that no changes have occurred, they could easily check this by\nmeasuring the different moulds and comparing whether the measurement on mould 1, 2 and 3 are all the same, which would\nvalidate the point they are trying to make.\n\nAuthor’s response to reviewer: The reviewer is correct that bioerosion can impact the measures collected here, which is why\nmarks where bioerosive processes impacted depth were not included in the analysis. The specimens cited above had\nmultiple modifications present, and those that were altered were excluded from our results.\nReviewer’s reply: have the authors clearly stated this in the text? It would be very easy to add a sentence that clarifies this. If\nthis is stated in the text, please can they indicate where, because I cannot find it.\n\nAuthor’s response to reviewer: We have added the information for the objective used in our analysis.\nReviewer’s reply – thank you for including this information\n\nAuthor’s response to reviewer: … For example, measuring a mark with a 20x lens would produce a point cloud that would\nexceed the memory capabilities of the S-Neox and would take hours instead of minutes to capture, if it were even possible…\nReviewer’s reply: I won’t go too much into this, but the idea that the point cloud exceeds the memory of the computer is\ninteresting... I, and other colleagues, study cut marks, tooth marks and percussion marks at 20x, and have no problem at all\nwith memory. Yes it takes longer, but science requires sacrifice, and while the work has been done on moulds, there is no\nissues with limited access to the samples which would imply that the authors could not have taken the time to scan all the\nmarks at 20x. Maybe older versions of the same microscope present storage issues that newer version of the S-Neox do not\npresent, if that is true then it might be necessary to state the year, version, or some indication about which edition of the\nmicroscope is being used, because I repeat, in our lab (University of Bordeaux), we have an S-Neox and have had no\nproblem studying cut marks at 20x magnification.\n\nAuthor’s response to reviewer: the quantitative method used effectively discriminates between the processed marks so\nremoving form is not problematic for our particular method and is necessary to capture the variables of interest.\nReviewer’s reply: I’ll have to take the author’s word for it, however they still have not stated anywhere which polynomial is\nused to remove the form from these scans. This is an incredibly simple statement to include in the methods – A polynomial of\nX was used to remove the form.\n\nAuthor’s response to reviewer: While there is some literature that agrees with [QDA should not be used for ordination] it\nmostly relates to LDA and geometric morphometrics and not QDA on variables related to size rather than shape.\nReviewer’s reply: There is a lot to unpack, even in this short paragraph. First a point of debate: “variables related to size\nrather than shape” - Size and shape can be, and are often, inherently related - unless you provide evidence to say otherwise,\nthis statement can be very problematic. Because the objective of the paper is not to look at this question, I will brush past this\nstatement, but I just would like to point out that care is needed here because this is not necessarily true, and in any other\npaper I would very strongly disagree with such a statement. Regarding the much more important point, I think the author’s\nhave missed the point with reference to the issues of using discriminant scores for ordination. I am not talking about\ngeometric moephometrics (GMM), I am talking about mathematics. All the literature I can find on QDA and the corresponding\nplots still relate them to the calculation of canonical variates. In fact, QDA and LDA are built on some fundamentally identical\nalgebraic calculations, that means they are completely related. Considering how QDA and LDA are heavily linked\nalgebraically, this would in turn imply that all issues related with LDA, which are related to CVA, are related in turn to QDA\nas well. I’m talking exclusively about algebra, not about GMM. The reason I mention GMM is simply because this is where\nthe majority of research into these algebraic properties lie. However, I would like to point out that authors such as Cardini and\nBookstein do not even use landmark based data in their papers about this topic, they use simple multivariate normal\ndistributions, and thus still prove my point – QDA should not be used for ordination. Likewise, the paper by Albrecht (that I\nprovide below) talks mostly about linear measurements, not landmarks… therefore it is directly relatable to the present study.\nIt is interesting that the scores in the plots are also called canonical something. From this perspective how can you be so\nsure the plot is different from what has been criticised by other authors? Likewise, no indications or citations are included in\nthe paper to clarify how calculations are performed, beyond the statement that JMP was used. JMP is a private licensed\nsoftware, which I do not have a licence for (in my opinion this in itself is an issue) so I cannot perform any tests to check for\nreplicability, or even evaluate the suitability and reliability of the algorithm in this specific context.\n\nIn sum; (1) if QDA is less sensitive to sample imabalance, i would like some references to support this, or simulations to\nprove this. (2) You still need to clarify if cross validation is done, and if it isn’t performed, then justification as to why must be\nprovided, because in a lot of circumstances the exclusion of CV would be considered completely problematic. (3) I am very\nhappy to see a PCA is now included, which would generally resolve a number of the issues I had with the QDA plot,\nhowever, there are a number of issues related to the PCA. (3a) The PCA provided is very very hard to read, and I cannot see\nwhere the fossil points fall. Have they even been included? Please improve the presentation of this figure so it is clear to\nread and see how the archaeological samples compare with the experimental ones. This is something you did with the QDA\nplot (you already know my opinion about how reliable this plot is), but have refrained to do for the PCA; (3b) why have the\nconfidence intervals been excluded for the PCA, but included for the (problematic) QDA? Please standardise the\npresentation of your results. (3c) Having red and orange points on the PCA makes the samples really difficult to distinguish.\nAlso, is it just me, or do the colours not match up between the QDA and the PCA? Or even within the PCA at all. In the key\nthere is blue, purple, green and purple, but in the PCA I cannot see blue, but instead see red, orange, green and purple.\nWhat is going on here?!\n\nIf the authors are still adamant that their use of QDA is appropriate, please provide justification. I provide justification here to\nhighlight why it is not appropriate (once again, before the authors claim that literature originating from GMM is not relevant\nhere, please read the methods sections of these papers to see that in most cases, the methods are not rooted in GMM\nanalyses at all);\nAlbrecht, G. H. (1992). Assessing the affinities of fossils using canonical variates and generalized distances. Journal of Human\nEvolution, 7(4), 49–69.\nBookstein, F. L. (2019). Pathologies of between-groups principal components analysis in geometric morphometrics.\nEvolutionary Biology, 46, 271–302.\nCardini, A., & Polly, P. D. (2020). Cross-validated between Group PCA scatterplots: a solution to spurious group separation?\nEvolutionary Biology, 47, 85–95.\nCardini, A., O’Higgins, P., & Rohlf, F. J. (2019). Seeing distinct groups where there are none: spurious patterns from\nbetween-group PCA. Evolutionary Biology, 46, 303–316.\nMitteroecker, P., & Bookstein, F. (2011). Linear discrimination, ordination, and the visualization of selection gradients in\nmodern morphometrics. Evolutionary Biology, 38, 100-114.\nRohlf, J. F. (2021). Why clusters and other patterns can seem to be found in analyses of high-dimensional data. Evolutionary\nBiology, 48, 1–16.\nIts incredibly important to point out that in each of these referenced papers, the authors may talk a lot about GMM, however\nthe data they are actually using to explore the algebraic and mathematical issues presented are simply multivariate\ndistributions. This means the data used for these statistical and theoretical experiments are completely analogous with the\npresent author’s data. In particular I recommend the paper by Cardini and Poly (2020) for proof on the importance of cross-\nvalidating their discriminant function plots. For the importance of cross-validation in classification I could recommend a large\nnumber of papers and books, however would simply state that the book by Hastie, Tibshirani and Friedman (2009) The\nElements of Statistical Learning is more than a good starting point.\n\nReply to the reviewer: We employ QDA here due to the unequal variance and size of the samples and we assess our model\nsetting prior probabilities proportional to the occurrence of each mark to account for the differences in sample size.\nReviewer’s reply: Ok then please clarify this in the text, and also state what prior probabilities were used. This can be easily\nincluded in a single sentence.\n\nReply to the reviewer: We have added an ROC curve as supplemental materials that reflects the requested information and\nshows the model is highly effective. Thanks for the suggestion! We will also report accuracy as it is a standard that more\nscientists are familiar with.\nReviewer’s reply: While the authors have already corrected the presentation of results to include more than just accuracy\n(ROC curves are a great addition), I have to simply state, for the purpose of scientific clarification, that I strongly disagree\nwith the term or concept that accuracy is a \"standard that more scientists are familiar with\" - I could not disagree more with the\nlogic that just because everyone else does it means that it is right. For years it was general consensus that the world was\nflat, this doesn’t mean that these individuals are the future of science. For years nobody thought anything was smaller than\nthe atom, scientific research proved them wrong too (if anyone is paying attention they might catch a Friends reference in\nhere). Science evolves and so should we, just because \"scientists\" (this should actually be\narchaeologists/palaeoanthropologist/palaeontologists) are more familiar with accuracy, that doesn’t mean we shouldn't\nbreak the mould and develop our analyses to adopt correct practices. Moreover, if you actually consult the literature\nregarding classification in data science, the authors would see that (mostly) non-archaeological scientists publish far more\nthan simply accuracy. Given this last point, I would argue that in data science reporting only accuracy is certainly not a\nstandard, which would make this statement incorrect on a number of different levels.\n\nJust to prove the point, and make sure that we go beyond the field of data science and into the realm of archaeology:\nimagine we have a new set of 100 marks; 80 cut, and 20 trampling. Imagine our algorithm claims that all marks are cut\nmarks. This would result in 80% accuracy, however, it would also be misleading considering how the algorithm is actually\nallocating the label cut mark to everything (i.e. hasn't learnt to classify anything at all). Now, if we calculate sensitivity and\nspecificity, at least one of these values would be 0, successfully helping us detect that our 80% accuracy is misleading, and\nthe algorithm is in-fact awful, and should not be used in real case studies (EVER!). Now imagine that instead of\narchaeologists we are doctors and we have potentially killed 80 patients. Just because our patients have been dead for 1.9\nmillion years (wow, what a date!), doesn't mean we should not be just as rigorous and strict with our diagnostics. I apologise\nfor rambling, but I find this point fundamental to clarify, especially if the authors find more great finds in the future and have\ndifferent reviewers who are unfamiliar with this concept.\n\nFinal statement - Once again, the reasoning for my lengthy review letter is to ensure clarity in my concerns, or open potential\ndebate for future reference, and I hope the authors find these comments constructive (if my tone is taken at some points to be\naggressive, I sincerely apologise, this is certainly not my intention). I strongly value and respect the work of all the authors,\nand think that these finds are incredible; I just want to ensure that there is no wiggle room or overlooked error that can be\nused against or to discredit this study. If I can find these little problems, then other (less amenable) researchers certainly can\nas well. I am incredibly grateful to be part of the review process, and look forward to seeing this paper published as soon as\npossible! Signed - Lloyd A. Courtenay", "reviewer_3": "I have carefully reviewed the authors' responses to my comments, the revised version of the manuscript, and the\nsupplementary information provided. I am grateful that they have considered my suggestions and that they are incorporated\ninto the manuscript.", "reviewer_4": "Hominin presence in Eurasia by at least 1.95 million years ago, by Curran et al.\nThe authors have convincingly addressed most of my comments, and even when objecting, the answers given are\nreasonable and make sense. After going through the detailed rebuttal, it seems to me that the evidence provided on the\nthree main critical aspects of work (i.e., human origin of the cut marks on bones; (ii) antiquity of the fossil material and (iii)\npalaeoenvironment reconstruction) are robust enough to support publication. I would especially like to thank the authors for\ndealing with the comments made by the reviewers in such a constructive way.\nThere are a couple of general comments I would like to make here after reading the revised manuscript, although they\nshould not preclude the publication of this very interesting work that undoubtedly deserves to be shared with the community.\nFirst, I must admit that I remain quite sceptical about the precision of the general chronostratigraphic framework of the ORV\nand Dacian basin. Yes, biochronology can give a rough chronological estimate of a given site, but this should be rather\nregarded as what it is, a rough estimation. Estimation because the inherent uncertainty typically associated to biochronology\nis usually minimized (or ignored). It may be related to (i) the taxonomic attribution of fossil specimens, which may vary\ndepending on the specialist involved, (ii) the validity/robustness of large scale biostratigraphy correlations across a country\nor even a continent, which may often be questionable, and (iii) and the existing true independent age control available at\nother sites used for comparison (i.e., based on methods like magnetostratigraphy and numerical dating method), which is\noften limited and has an uncertainty on its own. This is why I personally have difficulties to accept precise age constraints\nbased on biochronology only (e.g., here: ‘This formation extends from the base of the Pleistocene to as young as 1 Ma;\nBiochronological assessments indicate Grăunceanu is Late Villafranchian (~2.2-1.9 Ma) and is attributed to mammalian\nbiostratigraphic zones MN17/MmQ19-11). For example, in absence of numerical age results (and even though, the existing\nage errors are usually not negligible) or magnetostratigraphic evidence indicating the undisputed presence of Jaramillo, the\nminimum age constraint of 1 Ma proposed here is questionable. Let’s say that an Early Pleistocene chronology may roughly\nbe estimated for the Tetoiu Fm based on biostratigraphic evidence. Following a similar reasoning, there are many examples\nof localities that have been equally attributed to either the Middle, Late or Epivillrafranchian biochrons depending on the\nauthors. For example, I note that on the figure given in the rebuttal, the chronological window proposed for the site covers the\nlatest Middle Villafranchian and earliest Late Villafranchian, whereas the text mentions a much finer age window (‘Late\nVillafranchian (~2.2-1.9 Ma)’). Anyway, my point here is that I feel that the general chronostratigraphic context of the ORV\nand Dacian basin (and mainly the Tetoiu Fm) would deserve a dedicated section in Supplementary Material, as\nbiostratigraphy provides a crucial support to the U-Pb dating results. This would help the reader unfamiliar with the area\ngetting the bigger picture. In particular, it is unclear how the various sites mentioned in the text do relate (e.g., La Pietriș,\nValea Roșcăi, Fântâna lui Mitilan). It is mentioned that some of them are coeval, while other are younger/older, but no data\nare provided to support this. I understand that his information might be found in various other papers, but it would be\nbeneficial to the present work if a synthesis could be given in Supplementary Material. In particular the figure of the rebuttal\nshowing the biochronological evaluation of the Grăunceanu should probably be added. It does provide useful information.\nThe other localities form the ORV are that are mentioned in the text could be added to that figure too.\nFinally, I would like to add that, in my view, there is still an inherent general uncertainty on the context of the site and of the\nfossil assemblage. I do get the explanations of the authors on this matter, and I accept them, but the fact that the site is no\nlonger accessible, and that most of the information related to the excavation is gone, naturally raises some questions on the\nintegrity, extension and stratigraphy of the fossil assemblage. This is complicated by the fact that various successive\ncampaigns have been performed over more than a decade, and we ignore how they relate. Moreover, dating analyses were\nnot performed on the same bones that have the cut marks, or on the equid teeth used for palaeoenvironment reconstruction,\nso there is also an inherent uncertainty on whether all these remains are truly associated, i.e., come from the same\nexcavation area and the same stratigraphic unit. In that regard, would it be possible to at least know from which excavation\ncampaigns were taken all the specimens analysed here for cut marks, dating and paleoenvironement analyses? That said, I\nget the authors’ explanations and especially regarding the taphonomy suggesting a single depositional event, and no doubt\nthey did their best to narrow down this uncertainty, but this is still indirect evidence that cannot fully compensate for the\nabsence of a proper stratigraphy and spatial location of the remains.\nSpecific comments\n1. L.74-75: The chronostratigraphy of the Tetoiu Fm described from Andreescu et al. (2013) is not as robust and precise as\nsuggested here. There is no independent age control available, and if I am not wrong the age constraints are only based on\nbiostratigraphy. I’d rather see something like: ‘An Early Pleistocene age may be considered for this Fm based on\nbiostratigraphic evidence.’. I did not see evidence in Andreescu et al. (2013) for such a precise minimum age boundary of 1\nMa.\n2. L.104-106: Please indicate that this relative chronology is based on biostratigraphic evidence.\n3. Fig. 1: please check the spelling of ‘El-Kherba’.\n4. L.196-198: I understand these ages only refer to Grăunceanu. This should probably be clarified.\n5. L.199, ‘cluster tightly’: maybe a number (e.g., coefficient of variation, or a time range between the youngest and oldest\nresults) can be provided to support this statement.\n6. Fig. 3: please indicate in caption that age errors are 2σ.\n7. L. 231: VSMOW should probably be spelled out for its first occurrence in the text.\n8. Table 1 and Supplementary Table S7: The various ages reported there are quite confusing, maybe this should be clarified\nin caption. Sometimes an age range is given, with either 1 or 2 decimal digits. Not sure a 2 decimal digit precision is\nappropriate for most sites (Fântâna lui Mitilan even shows 3 digits). It is also unclear whether the age range is taken from the\nerrors on the numerical ages, from the combination of various ages, or simply correspond to various levels from the same\nsite.\n- Not sure to understand how the age range of 2.3-2.0 Ma for Grăunceanu was obtained since the various U-Pb ages range\nfrom 1.87 to 2.01 Ma\n- Ain Hanech and El Kherba (Algeria) are usually regarded as coeval. Please see Duval et al. (2023) for the most recent age\nestimation.\n- Ain Boucherit (Algeria) is made of two levels, the lower level (Lw) dated to ~2.4 Ma and the upper level (Up) to ~1.9 Ma.\nMaybe these 2 different levels (usually referred as Ain Boucherit Lw and Ain Boucherit Up) should be differentiated,\notherwise the age range given may suggest a significant uncertainty on the age of the site, which is not the case.\n- Pirro Nord: please note the most recent dating study on that site (Duval et al., 2024). Although a biochronological age of\n1.3-1.6 Ma cannot be excluded, there is strong evidence suggesting that the site may be much younger.\n- Sima del Elefante : maybe indicate the layers, e.g., Sima del Elefante (TE7-9). Methods also include BS.\n- Fuente Nueva-3: a rounded age of ‘1.2 Ma’ would probably be more appropriate given the associated uncertainty.\n- Barranco Leon: Methods also include MS and BS, in addition to ESR.\n9. L.324-330: for a matter of consistency, please use the same number of decimal digits for all percentages.\n10. L.326: Do the 3.32% include both Lower and Upper levels at Ain Boucherit? Please clarify.\n11. L.335: maybe you can refer to Supplementary Table S8 here.\n12. Ref 71 in SM: please note the name of the journal is ‘L'Anthropologie’.\n13. Supplementary Table S8: please indicate in caption the meaning of ‘Full Taph?’. Note that the age reported for\nGrăunceanu is not consistent with that in Table S7. Same comment applies for Dmanisi.\nReferences\nDuval, M., Arnold, L.J., Bahain, J.-J., Parés, J.M., Demuro, M., Falguères, C., Shao, Q., Voinchet, P., Arnaud, J., Berto, C.,\nBerruti, G.L.F., Daffara, S., Sala, B., Arzarello, M. (2024). Dating the earliest evidence of human presence in western Europe:\nNew results from Pirro Nord (Italy). Quaternary Geochronology 101519. 10.1016/j.quageo.2024.101519.\nDuval M., Sahnouni M., Parés J.M., Zhao J.-x., Grün, R., Abdessadok S., Pérez-González A., Derradji A., Harichane Z.,\nMazouni N., Boulaghraief K., Chelli Cheheb R., van der Made J. (2023). On the age of Ain Hanech Oldowan locality\n(Algeria): First numerical dating results. Journal of Human Evolution 180, 103371.\nhttps://doi.org/10.1016/j.jhevol.2023.103371.\nSahnouni M., Parés J.M., Duval M., Cáceres I., Harichane Z., Van der Made J., Pérez-González A., Abdessadok S., Kandi\nN., Derradji A., Medig M., Boulaghraif K., Semaw S. (2018). 1.9-2.4 million-year-old artefacts and stone tool cutmarked\nbones from Ain Boucherit, Algeria. Science 362 (6420), pp. 1297-1301. DOI: 10.1126/science.aau0008."}, "author_response": "We appreciate the reviewer’s continued support of our manuscript and their thoughtful\ncomments on the statistical methods. It is clear from their comments that they feel strongly\nabout these methods, and we have done our best to address their concerns in as thorough a\nmanner as possible. We agree that some of the confusion here may have been related to\ndifferences in terminology, and we agree that a robust cross-validation of the model is critical\nfor ensuring proper interpretation of the results.\nTo that end, we have revised the text to make our analyses as clear as possible, and now\ninclude an additional author who performed updated analyses. As per the reviewer’s\nsuggestion we now present both leave-one-out and 10-fold cross-validation of the QDA\nresults. We have added a new Supplementary Table 3 that provides the confusion matrices,\nand have added the R code to run these analyses in the supplementary information (though\nnote this requires the raw experimental data, which are not yet published). Both cross-\nvalidation methods return very similar results (81.4% and 82%, respectively), demonstrating\nthe effectiveness of the QDA. In the course of these updates we noted a small error (n=3) in the\ncounts of bone surface modifications in the experimental dataset and have subsequently\nupdated Supplementary Figures 22 and 23, Supplementary Table 2, and all posterior\nprobabilities listed in the text. These changes are all very minor and do not change the\noutcomes of our interpretations; in fact, they strengthen them.\nAll changes made have either been tracked in the text or sections with major changes are\nhighlighted in yellow. Again, we appreciate the reviewer’s attention to our statistics and\nsupport of our work, and we hope these updates satisfy their concerns regarding the\nquantitative assessments of the bone surface modifications we present here."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors have finally acknowledged and made modifications in light of my concerns, which I am grateful for. I want to\nemphasize that my intention was never to be overly critical or harsh. Rather, I believe it is fundamental to evaluate\nquantitative analyses rigorously. My recommendations were aimed at reinforcing the contents of the paper, and my approach\nwas intended to be constructive, not deconstructive. As I have always stated, I had no issue referring to the marks as cut\nmarks. However, statistical analyses require specific criteria to be met, and I was simply trying to highlight those\nrequirements.\nI would like to thank the authors for all their hard work and patience. This is truly a fascinating and important study.\nKind regards\nLloyd"}, "author_response": ""}]

Decision letter