Hominins on Sulawesi during the Early Pleistocene

Palaeomagnetic analysis

Oriented block samples of approximately 10 × 10 × 10 cm³ each were carved with a sharp non-magnetic knife from fresh excavation baulks (samples Sp1 and Sp2) or from field exposures for palaeomagnetic analysis (Sp3–Sp14). Strike and dip were measured using a Brunton Compass corrected for local declination. Owing to the scarcity of suitable consolidated fine-grained layers, some samples were obtained from outcrops located up to 430 m from the measured section along the main Cabenge–Pampanua road (Extended Data Fig. 1). This stratigraphic section was originally recorded in 1991, when the outcrops along the road were still freshly exposed by construction works¹⁹. Five samples (Sp1–Sp5) were taken west of the Bulu Cepo Fault in subhorizontal to slightly west-dipping bedding, including two samples from the Calio excavation (Sp2 and Sp3), one sample from an outcrop along an erosional gully 74 m northwest of the Calio excavation (Sp1), one sample from the Beru (modern Indonesian spelling Berru) test-pit at 328 m eastsoutheast of the Calio Excavation (Sp4) and the fifth sample (Sp5) from an outcrop on the opposite side of the road 17 m south of the Beru test-pit. Samples Sp6 to Sp14 were obtained from steeply west-dipping layers east of the Bulu Cepo Fault. For those block samples taken at a distance from the road section, the approximate stratigraphic position relative to the section recorded in 1991 was extrapolated using the measured strike and dip of 258°/60° (Extended Data Fig. 1d).

The carved samples were tightly packed in plastic and shipped to Bandung, Java, where they were further analysed at the Palaeomagnetism Laboratory of the Centre for Geological Survey (CGS). One sample (Sp4) arrived in a deteriorated condition and could not be analysed further. Of the remaining block samples, three to four 1 × 1 × 1 cm³ subsamples (denoted specimens and labelled sequentially) were carved from each, or, in cases in which the samples were not well consolidated, plastic cubes of the same dimensions were pressed into the sediments. Subsamples were stored in magnetically shielded containers to minimize viscous remanent acquisition of any modern magnetic field components. Care was taken to keep samples moist until palaeomagnetic measurements were made at the CGS Palaeomagnetic Laboratory.

Of all the palaeomagnetic block samples collected, three sequentially numbered subsamples (specimens) were then step-wise demagnetized with a Schonsted USA GSD-1 series alternating-field demagnetizer (AFD), at intervals of 2.5–5 mT, up to peak fields of 60–100 mT. If a fourth specimen was available, it was also step-wise demagnetized by thermal demagnetization (TD) using a magnetic measurement thermal demagnetizer MM80 series, at intervals of 100–25 °C, with a temperature up to 800 °C. Magnetic remanence directions and intensities were determined using a dual spinnermagnetometer from AGICO (model JR-6A). Most of the samples were influenced by a secondary magnetization; however, this could be easily removed. Characteristic remanent magnetizations (ChRMs) were isolated throughout stepwise AFDs at 20–40 mT in level, above which most samples were completely demagnetized (Extended Data Fig. 7). Meanwhile, ChRMs were isolated through stepwise TDs at 200–300 °C to 550–600 °C, above which most samples were completely demagnetized (Extended Data Table 1 and Extended Data Fig. 8). The temperature range for ChRM indicates that magnetite is the dominant carrier of remanence. The ChRM directions obtained by TD are consistent with those obtained by AFD (Extended Data Table 1).

The ChRM is either trending to the origin of the orthogonal vector projections or is defined as the mean of vectors which are stable. Fisher’s statistic was used to establish the mean direction of primary magnetism (Extended Data Table 2). The primary magnetic direction was determined from at least the last four to five plot points using principal component analysis³¹, with the mean maximum angular deviation set at <15° using PuffinPlot software³² and IAPD 2000 (ref. ³³). Magnetic polarities are defined from the virtual geomagnetic pole (VGP) where a normal polarity represents positive (Northern Hemisphere) VGP latitudes and a reverse polarity represents negative (Southern Hemisphere) VGP latitudes. Intermediate or transitional polarity fields are defined as those with VGP latitudes between 45° N and 45° S. Each VGP is calculated using the unit vector mean of the three to four directions observed at each sampling level.

Extended Data Fig. 9 shows the VGP plot of the 13 samples from the Calio excavation site and the samples taken further east. Twelve of the VGPs cluster in the Southern Hemisphere around the southern part of South Pacific Ocean and Antarctica, and all exhibit a reverse polarity. One sample (Sp9) has an intermediate or transitional VGP latitude of 38.8° S, while another sample (Sp7) has a VGP latitude on the Northern Hemisphere (normal polarity).

The four samples associated with the subhorizontal fluvial beds of the Calio excavation site and Beru-Bulu Carulle test-pit west of the Bulu Cepo Fault all have a reverse polarity. This suggests that the sampled horizons have a minimum age of 773 ka, which corresponds to the Brunhes–Matuyama boundary²⁵. On the basis of the US–ESR dating of the two Celebochoerus teeth, the approximately 12-m-thick sequence sampled at the excavation site and surrounding locales is more likely to pertain to the reverse magnetic interval below the Jaramillo Subchron (1,070–990 ka)²⁵.

Three samples taken from the steeply west-dipping flank exposed at Ciangkange (Sp6, Sp10–11) also have a reverse polarity, while one sample (Sp9) shows an intermediate polarity (Extended Data Fig. 5). This interval is probably still within the Matuyama Chron. The normal polarity of sample Sp7 further down the sequence could be correlated to the Olduvai Subchron dated to 1.948–1.787 Ma (ref. ²⁵), but this would need to be confirmed by a higher sampling density in the stratigraphy.

U-series dating

U-series data were obtained by first sectioning a part of each of the two Celebochoerus fossil teeth using a high-precision diamond blade (300 μm width) to expose both dental tissues. The exposed surface of each sample was then polished to >10 μm smoothness to provide a clean ablation surface. U-series measurements were undertaken by laser-ablation multi-collector-inductively coupled plasma mass spectrometry (MC-ICPMS) at the Geoarchaeology and Archaeometry Research Group (GARG) Biomics facility at Southern Cross University. Laser ablation was performed using a New Wave Research 193 nm Arf excimer laser, equipped with a TV2 cell. Thorium (²³⁰Th, ²³²Th) and uranium (²³⁴U, ²³⁵U, ²³⁸U) isotopes were measured on the Thermo Neptune XT MC-ICPMS system mounted with jet sample and x-skimmer cones. All five isotopes were collected in static mode, with both ²³⁴U and ²³⁰Th collected in the ion counter and CDD respectively. Helium flow rate (850 ml min⁻¹), nitrogen (6 ml min⁻¹) and ICP-MS parameters were tuned with NIST610 element standard to derive a ²³²Th/²³⁸U ratio of greater than 0.85 and therefore minimizing differences in fractionation between Th and U. Tuning was achieved with a fluence of about 6.5 J cm⁻², pulse rate of 100 Hz, spot size of 120 μm and scan speed of 5 μm s⁻¹, yielding 4.40 V of ²³⁸U and 3.85 V of ²³²Th on NIST610.

The Celebochoerus maxillary teeth were ablated using rasters of ~5 min each (twice ~800 μm long). Before and after each sample, NIST612, MK10 and MK16 standards were measured³⁴, as well as a fossil rhinoceros tooth with known isotopic ratios³⁵. The ²³⁴U/²³⁸U and ²³⁰Th/²³⁸U isotopic ratios were corrected for elemental fractionation and Faraday cup/SEM yield by comparison with MK10 coral for which ratios were previously characterized internally by solution analysis³⁴. Detrital-corrected ²³⁰Th-U ages were calculated for each analysis using IsoPlotR³⁶ and UThwigl³⁷ with an assumed detrital (²³⁰Th/²³²Th) activity ratio of 0.8 ± 0.8 (ref.³⁸). Concentrations of U and Th were determined using NIST612 glass as a calibration standard. Background subtraction, concentration quantification and ratio corrections were performed using Iolite software³⁹. The corrected ²³⁴U/²³⁸U and ²³⁰Th/²³⁸U isotope ratios for the secondary standard (MK16 coral) within error of the value were determined by solution analysis. The rhinoceros tooth was used as a control for matrix effect.

Little to no detrital thorium was measured within the dental tissues, except at the surface and along cracks containing sediment intrusion (Extended Data Table 3). For the measurements, we avoided macroscopically visible cracks and obvious altered zones. The distribution of uranium across the dentine and the enamel is characteristic of the uranium diffusion pattern described in ref. ⁴⁰. Uranium diffused from the pulp cavity across the dentine and accumulated at the enamel–dentine junction before diffusing slowly into the enamel crystal structure. This translates into a large variability of uranium content at the interface of the two tissues (Extended Data Table 3). In summary, the Celebochoerus teeth from Calio exhibit varying uranium distribution patterns (Extended Data Table 3), including heterogeneous age distributions, as well as uranium accumulation in altered zones and along the enamel–dentine junction. Such values across dental tissues are characteristic of the complex diffusion history of ancient fossil teeth with several diffusion episodes, localized uranium hotspots and leaching zones, especially close to the pulp cavity and along cracks, as well as other alterations in the dental tissues⁴⁰.

US–ESR dating

The ESR dating was carried out on the Celebochoerus fossil teeth (left P⁴ and left M³). The outer layer of enamel was extracted mechanically from dentine with a rotary saw and cleaned carefully of any dentine residue (outer ~100 µm ± 10% on each side). The P⁴ tooth enamel was then ground into powder and sieved to 90–180 μm before being split into 10 aliquots (about 95–100 mg each). While one aliquot was kept as a natural dose, the other nine aliquots received gamma irradiation steps of 50, 100, 250, 600, 1,200, 2,400, 4,000, 8,000 and 15,000 Gy at the Australian Nuclear Science and Technology Organisation. The aliquots were irradiated by a calibrated Gammacell 220 (⁶⁰Co) gamma source with a dose rate of ~23.8 Gy min⁻¹. The M³ enamel was kept as a fragment, and this fragment was mounted into a parafilm mould within a Teflon sample holder to record the angular dependency in the ESR response^41,42,43,44. Irradiation of the fragment was performed with the Freiberg X-ray irradiation chamber for fragments, which contains a Varian VF50 X-ray gun at a voltage of 40 kV and 0.5 mA current on the fragment exposed to X-rays without shielding (apart from a 200 µm aluminium foil layer)^45,46. Each fragment was irradiated after exponentially increasing irradiation times (at 90 s, 380 s, 900 s, 1,800 s, 3,600 s, 7,200 s, 14,400 s). For each irradiation step, the energy output of the X-ray gun was recorded at the beginning and end and then averaged, enabling us to correct for the dose rate received. For each irradiation step, the fragment was measured over 180° in x, y and z configurations with a 30° step^41,43.

ESR intensities of tooth enamel powder were calculated using calibrated quartz tubes and measured using a Frieberg MS5000 ESR spectrometer at the GARG facility. Each tube was rotated three times with a 120° increment. All measurements (powder and fragment) were performed with the following conditions: 2 mW microwave power, 100 kHz modulation frequency with 1,024 points resolution, 0.1 mT modulation amplitude, 45 s conversion time, 12 mT sweep width and 21 s sweep time. The ESR intensities and dose–response curves followed previously reported recommendations²⁹. The D_E was determined using the McDoseE 2.0 program⁴⁷. Both single saturated exponential (SSE) function and double saturated exponential (DSE) function were applied to obtain D_E values, and both functions provided statistically undistinguishable results (P⁴: DSE, 2,275 ± 124 Gy; and SSE, 2,267 ± 99; M³: DSE, 2,140 ± 108 Gy; and SSE, 2,116 ± 118). Yet, it is our experience that SSE provides more systematic and reliable results than DSE, as long as recommendations from ref. ²⁹ are followed.

ESR intensities were extracted from T1–B2 peak-to-peak amplitudes on the merged ESR signal or the powder spectra⁴¹. Isotropic and baseline corrections were applied uniformly across the measured spectra for fragment measurements⁴². The amount of unstable non-orientated CO₂ radicals (NOCORs) was estimated on the fragment irradiated by X-ray to be negligible (~1%) and, therefore, could not be applied to the powder sample (M³) using the protocol described previously⁴¹. The results of the gamma rays irradiated powder are consequently uncorrected for NOCORs. The ESR dose–response curves for the fragment were obtained using merged ESR intensities of all orientations and associated s.d. values from the repeated measurements over one orientation only. For powder spectra, associated s.d. values were obtained using repeated measurements after shaking the tube at each irradiation step.

Sediment elemental concentrations, external beta and gamma dose rate contributions, and water content (measured) are shown in Extended Data Table 4. The external beta dose rates have been extrapolated from the U, Th and K contents measured on the associated sediment (~50 g). The sediments were collected in the immediate vicinity of the sample. The sediment was crushed, homogenized and separated into three aliquots before being digested in a 1:3 nitric/HCl acid solution (APHA 3125 ICPMS). The external gamma dose rates were determined using the associated sediment and assuming a 4π geometry. The cosmic dose rate was estimated based on consideration of the altitude, geomagnetic latitude and density (2 g cm⁻³) of sediment overburden⁴⁸. However, a large part of the ancient sedimentary deposit (probably several metres in thickness) at the Calio locality has clearly eroded away over time, leaving only a 30 cm depth for the current burial sample.

The internal dose rate was calculated from the U-series measurements obtained on each enamel fragment and surrounding dentine. Both enamel and dentine were measured using rasters (see the U-series dating description above for parameters). Concentration and isotopic ratios were measured directly on the fragment used, and all measurements were averaged to obtain a single value for the entire tooth. Post-Rn equilibrium was assumed. Ages were calculated using dose rate conversion factors described previously⁴⁹. Moreover, ESR ages were calculated using both the US-ESR MATLAB program⁵⁰ and the DATA program⁵¹ to compare US–ESR results, and with DATA to obtained closed-system calculation (CSUS-ESR age estimation).

Stone artefact analysis

The stone artefacts excavated from Calio were analysed using qualitative and quantitative methods published elsewhere^52,53,54. Artefact condition was assessed according to criteria outlined in Shea⁵⁵, with taphonomic abrasion and edge-damaged considered separately. Artefact material is classified as chert in all cases pending closer assessment of local stone sources.

Technological length refers to the measurement from the point of force application (PFA) to the distal end of the flake along the percussion axis. Technological width is the maximum width at right angles to technological length. Thickness refers to maximum thickness at right angles to the plane defined by technological length and width. Maximum length and width are the maximum dimensions of the flake, along any axis. Thickness was also measured at 25%, 50% and 75% of the distance from the PFA to the distal end, as a proxy for morphological uniformity.

Platforms describe the nature of the surface marked by the PFA. Platform depth is the distance from the PFA to the dorsal platform edge, as measured along the plane defined by the percussion axis. Platform width measures the maximum width of the flake’s platform⁵⁶. The number of flake scars on the platform surface was counted.

Exterior platform angles were measured using a lockable sliding bevel gauge combined with a digital goniometer. One arm of the sliding bevel was placed along the dorsal flake surface and the second arm was placed along the platform surface, from the dorsal platform edge to the PFA. The bevel was locked and the angle between the gauge’s arms was measured using the digital goniometer. The interior platform angle was measured directly using the digital goniometer. One arm was placed along the platform surface from the dorsal platform edge to the PFA, and the other arm was placed along the surface within 1–2 mm below the PFA, along the side of the hertzian cone at the start of the bulb of percussion, as described previously⁵⁷.

The sizes of retouching scars measuring >2 mm in maximum dimension were recorded. Length was measured from the existing platform edge to the distal edge, along the percussion axis. Width was measured at right angles to length. Cortex proportion was assessed visually as described previously⁵⁸. Cortex on platforms was included in the assessment of cortex coverage. The location of cortex was recorded as described previously⁵⁹.

Flake types were recorded as described previously⁶⁰, where type IV refers to flakes with non-cortical platforms and 100% dorsal cortex coverage; type V refers to flakes with non-cortical platforms and partial dorsal cortex coverage; and type VI refers to flakes with non-cortical platforms and no dorsal cortex (types I–III were not present in the Calio assemblage). Dorsal scar orientations followed the sector approach described previously⁶¹, excluding the redirecting flake.

Annotated 3D photogrammetric models were made of the Calio artefacts according to protocols described previously⁶². The photographs were taken using a Nikon D7200 camera with AF-S DX Nikkor 18–140 mm lens and processed with Agisoft Metashape v.2.0.1.16069. The 3D models can be accessed online (https://une.pedestal3d.com/r/ZdPcP-kXvc/).

The Calio artefact data are presented in Supplementary Table 1. Incomplete linear dimensions are recorded as negative values and highlighted in red font.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.