Hunter-gatherer sea voyages extended to remotest Mediterranean islands

Overview

Our multidisciplinary study combines archaeobotany (Supplementary Table 1), chronological modelling (Supplementary Tables 2–6), isotopes (Supplementary Table 7), contextual research and broad regional chronological modelling (Supplementary Tables 8 and 9) with anthracology and phytoliths (Supplementary Tables 10–13), lithics (Supplementary Table 14) and the study of faunal remains (Supplementary Tables 15 and 16). The methods used are described below, with further contextual information in the Supplementary Information and Extended Data Figs. 1–11.

Excavation and sedimentology

Here, we describe the excavation of a 5 × 5-m trench, designated Trench 4, at Latnija between 2021 and 2023, expanding on a 1 ×1-m test trench excavated in 2019. We set up an alphanumeric grid system in the doline to label each individual 1 ×1-m square, aligned in orientation with the 2019 test trench and with the nearby cave wall, with letters running on a SW–NE axis and numbers increasing on a NW–SE axis. The 2019 test excavation targeted square M2, with the expanded Trench 4 spanning squares J–N and 2–6, located at the northern edge of the doline spanning the dripline (Fig. 1).

Excavation was performed using a single-context recording methodology to resolve between discrete sediment units, with arbitrary subdivisions within a single deposit as 5–10-cm spits where necessary to aid control of find recovery and sediment sampling. Features of post-depositional disturbance, such as animal burrows, were readily differentiated from undisturbed sediments owing to their mixed character and friable texture and the presence of sediment voids, and were excavated in their entirety and excluded from our analyses. Finer-scale post-depositional disturbance occurs as limited fine rooting and is restricted to the uppermost deposits. The natural deposition of clasts from the shelter wall presents an alternate form of potential post-depositional disturbance that might have led to localized soft-sediment deformation. The three-dimensional position of all artefacts larger than 20 mm, bones larger than 20 mm and charcoal, and the geometry of excavation context boundaries, were recorded using a total station. Bulk sediment sampling retained a minimum of 60 l per context (predominantly in the uppermost deposits) up to 100% sampling of sediments, which were processed by bucket flotation using 250-µm mesh for macrobotanical recovery, followed by wet sieving through 5-mm screens for artefact recovery; sediments that were not retained for flotation and wet sieving were dry sieved through 5-mm screens. Additional sediment samples were recovered from each context for ancillary analyses.

So far, we have identified 309 discrete sedimentary contexts, reaching a maximum depth of 1.48 m from the surface. We have grouped contexts into six phases (Phases I to VI) on the basis of major changes in sediment colour, texture, composition and structure, alongside patterns evident in material culture. The stratigraphic matrix for the Mesolithic Horizon and immediately underlying deposits is presented in Extended Data Fig. 3.

Micromorphological samples for thin-section production were collected by cutting in situ, orientated blocks of sediment into Kubiena tins (90 mm × 70 mm × 50 mm). The location of Kubiena samples was dictated by the architecture of the sediment sequence and representative sediment deposits, and the contacts between deposits were targeted. The laboratory samples were air-dried for two weeks and placed in labelled plastic pots. The samples were immersed in a mixture of clear casting resin (four parts) to acetone (one part). To accelerate curing, a catalyst of methylethylketone peroxide was added (3 ml catalyst to 2,000 ml resin). The resin mixture was poured around the side of the sample to allow the larger pore spaces to be filled from adhesion and cohesion, and then completely immersed in the resin. The samples were impregnated under a stepped-vacuum regime to a maximum vacuum pressure of −25 in Hg for eight hours. The samples were left to cure for around six weeks until the resin was hardened, followed by a final cure at 65 °C for 15 h. The blocks were removed from the sample frame, split along their long axes and one surface polished on fixed diamond abrasives with successively finer grades (70 µm, 45 µm and 20 µm). The polished sample was stuck to a labelled slide using an epoxy resin that cures overnight. The slide and sample were cut down to around 1 mm and then excess sample was removed using a Jones and Shipman surface grinder. The sample was hand-polished to finish off the surface before coverslipping the sample again by bonding with an epoxy resin. Analysis of the thin sections was performed on a Leica M205C petrological stereo zoom microscope and image capture was done using the Image Pro-Express software.

Archaeobotanical methods

Studies of plants in the Mesolithic Horizon at Latnija were performed in the form of pollen analyses, anthracology, hearth phytolith analyses and macrobotanical identifications from remains recovered through flotation. These analyses were performed to reconstruct the vegetation of the site, determine whether any domesticated plants were present, investigate the use of different fuels at the site and unravel mineral composition to identify combustion structures.

For pollen analysis, we collected sedimentary samples to perform palynological analyses focused on the reconstruction of past vegetation at and near Latnija. Sampling was performed in Phase V contexts (034) and (048), both of which are characterized by the presence of thick ash and combustion residue deposits. This approach was adopted to correlate the palaeobotanical remains preserved in the sediment with human activities during Phase V, which is characterized by the oldest Mesolithic.

Samples were treated following pollen concentration techniques⁵². This included sediment deflocculation with sodium pyrophosphate, Lycopodium tablets with known content to calculate palynomorph concentration values⁵³ and 7-µm nylon sieve to discard clay-sized particles. Carbonates were removed with 10% HCl and concentrated at 2,500 rpm for three minutes. Heavy liquid separation using sodium metatungstate with a specific gravity of 2.0 and centrifugation at 1,500 rpm for 20 min was done to separate organic and mineral fractions. After recovering the upper supernatant fraction, this step was repeated to increase the concentration. The remaining fraction was treated with cold 40% HF for one night to eliminate remaining silicates. The residue was washed in 98% ethanol, glycerol was added and the remaining ethanol was evaporated. The solution was kept in glycerol, mounted on slides and identified at 400× magnification under a light-transmitted microscope by referring to established literature^54,55. Pollen counts were done up to 250 identifiable grains. A pollen diagram (Extended Data Fig. 4a) indicating values for each taxon as percentages of the total pollen sum was plotted with the help of C2 software⁵⁶.

For anthracological analyses, bucket flotation was used to recover charcoal and other carbonized archaeobotanical remains from the sediments, all of which were collected. Charcoal was also handpicked to provide a larger number to select for dating purposes and anthracology.

A total of 165 charcoal fragments were observed under reflected light microscopy (Motic PANTHERA) with dark and bright fields and ×50, ×100, ×200 and ×500 magnifications. Images were taken with an environmental scanning electron microscope (FEI Quanta 600) coating charcoal with gold. Each charcoal piece was manually fragmented into the three wood anatomy sections (transverse, tangential and radial). Observing the three anatomy sections allowed us to identify taxonomic characters. Different wood anatomy atlases and a comparative collection at the Catalan Institute of Human Paleoecology and Social Evolution were used to support the identifications^57,58. The assemblage is characterized by a number of indeterminable fragments related to wood anatomy alterations (cracks and vitrification) and/or size of the fragments.

To study the pyroarchaeological record of the Mesolithic Horizon at Latnija, we combined the study of phytoliths and the mineralogical composition of sediments by FTIR. We analysed 24 samples that were collected during the 2022 fieldwork from a large combustion structure identified in square N2 at the base of Phase V (Fig. 2). Sampling was performed on the basis of visual identification of the internal structure of the hearth, distinguishing between samples coming from the possible combustion residue (n = 10), samples coming from the thermal impact zone (n = 8) and control samples from below the hearth (n = 6).

Phytoliths were extracted following the fast extraction method⁵⁹. Phytolith quantification and identification was done using a Zeiss Axioscope transmitted light microscope at ×200 and ×400 magnifications. Phytolith morphological identification followed the standard literature and modern plant reference collections^33,60,61,62. We followed the terminology of the International Code for Phytolith Nomenclature (ICPN 2.0) for phytolith descriptions⁶³.

The mineral composition of the samples was identified using a Jasco FT/IR-6700 spectrometer. Infrared spectra were collected in the 4,000–400 cm⁻¹ wavelength range at a resolution of 4 cm⁻¹ using the conventional KBr pellets method. The spectra were interpreted using the position of the main peaks described on reference collections⁶⁴. Thermally altered clay was identified on the basis of specific absorption peaks in the clay spectrum⁶⁵, and the presence of anthropogenic or geogenic calcite was determined following previous studies^66,67.

The archaeobotanical samples from Latnija’s Mesolithic Horizon were recovered from the 2021 and 2022 excavation seasons. Although we engaged in a 100% sediment collection strategy, after flotation, not all samples from these phases contained plant macrofossils. The assemblage suitable for study consists of 28 samples in total—19 from the 2021 field season and 9 from the 2022 season. Each sample was processed in the field using a basic bucket flotation method, as described previously^68,69. The samples were then sent to the Max Planck Institute of Geoanthropology in Jena, Germany, for analysis. Once in the laboratory, samples were passed through nested U.S. Geological sieves to ease sorting. Material smaller than 0.50 mm was not sorted. Carbonized wood fragments larger than 2 mm were counted, although wood identification was done as a separate analysis and is reported above. Seeds and seed fragments were separated from all sieved contexts, and charred seeds were systematically collected. The identified taxa are presented in Supplementary Table 1.

Radiocarbon dating methods

Except for the bone samples, radiocarbon dating was performed at the Curt-Engelhorn-Centre Archaeometry (CEZA) in Mannheim, Germany. Samples included charcoal, seeds and marine shells. Bone samples were analysed at the University of Georgia Centre for Applied Isotope Studies (CAIS). We used a multistep chronological study to clearly constrain the Mesolithic Horizon at Latnija. First, we constructed a chronological framework for the site, which involved 31 charcoal samples and the one bone (Supplementary Tables 2 and 3). Charcoal samples were selected from contexts directly underlying the Mesolithic Horizon to help constrain the onset of Mesolithic occupation, excluding samples from burrows that appear at the interface of major divides in sediment depositional processes (Phases VI–V) (see Extended Data Fig. 5 for illustrated sample locations). In addition, charcoal samples were selected from contexts throughout the Mesolithic Horizon (Phases V–III), including direct sampling from hearths that appear at the base of Phase V (Fig. 2, Supplementary Information 2 and Extended Data Fig. 6). The model was divided into the major phases recorded during the excavation (Supplementary Information 2).

To obtain independent verification of the integrity of the age model, we also targeted marine gastropods (P. turbinatus in particular) because they formed clear in situ tip lines identified in Phase III. Forty-nine samples of P. turbinatus were dated for this purpose. The number of samples was chosen to reflect the fact that: (i) marine calibration is more complex than terrestrial calibration, thus a larger sample size was required to account for the natural spread in the data; and (ii) these shells are a direct measure of human presence, because they have been imported to the site by people.

Charcoal samples were prepared using a standard ABA pretreatment. This covers an acid step with diluted hydrochloric acid to remove calcite and lime attached to the sample. A base step with diluted sodium hydroxite follows to remove soluble humic acids. As the base attracts fresh CO₂, another acid step finalizes the pretreatment and removes any modern contamination. The samples are then combusted in an elemental analyser (MicroCube, Elementar) and the CO₂ is collected and graphitized to elemental carbon. The carbon is pressed into a target and measured in a MICADAS mass spectrometer⁷⁰.

The shell samples only undergo a treatment with diluted acid to remove adjacent carbon contamination from limestone or calcite. For shell samples, the CO₂ is extracted using phosphoric acid in an autosampler before graphitization, and measurements are the same as for the charcoal samples described in a previous study⁷¹.

The bone sample was cleaned by wire brush and washed using an ultrasonic bath. After cleaning, the sample was then reacted under vacuum with 1 M HCl to dissolve the bone mineral and release CO₂ from bioapatite. The residue was filtered, rinsed with deionized water and, under slightly acid conditions (pH 3), heated at 80 °C for six hours to dissolve collagen and leave humic substances in the precipitate. The collagen solution was then filtered to isolate pure collagen and dried out. The dried collagen was combusted at 575 °C in evacuated and sealed Pyrex ampoules in the present CuO. The resulting CO₂ was cryogenically purified from the other reaction products and catalytically converted to graphite. Graphite ¹⁴C/¹³C ratios were measured using the CAIS 0.5 MeV accelerator mass spectrometer. The sample ratios were compared with the ratio measured from the Oxalic Acid I (NBS SRM 4990). The uncalibrated dates were then given in radiocarbon years before 1950 (years bp), using the ¹⁴C half-life of 5,568 years. The error is quoted as one standard deviation and reflects both statistical and experimental errors. The date has been corrected for isotope fractionation. As with other terrestrial radiocarbon dates in this study, calibration was performed with OxCal 4.4 using IntCal20 and as part of a phase model for the site. Modelled and unmodelled calibrated dates and model diagnostics are presented in Supplementary Tables 2 and 3.

To correct for the MRE, we compared the ages of P. turbinatus shells with the ages of charcoal from the same stratigraphic contexts. The reservoir effect ΔR was modelled in OxCal 4.4. with the latest datasets of IntCal20 for the charcoal samples and Marine20 for the shells. It was modelled using a phase model and choosing a wide restriction for ΔR. Samples marked by OxCal as outliers are presented in the table but are not included in the next modelling step if the model cannot deal with them leading to an A of less than 60%. These outliers might reflect processes such as bioturbation. The results of the MRE calculations are shown in Supplementary Table 4, and the corrected dates for each P. turbinatus age are shown in Supplementary Table 5 (see also Supplementary Information 4).

Radiocarbon modelling methods

Models involving radiocarbon dates were used to address the key question of whether there is evidence of occupation in the Latnija cave excavation sequence that securely relates to human activity predating the available evidence for Neolithic habitation elsewhere on Malta and in the surrounding Mediterranean archaeological record. This was done by: (1) establishing the age of the Mesolithic deposits at Latnija; (2) determining when the wider regional Mesolithic-to-Neolithic transition is most likely to have occurred; and (3) determining whether there is evidence for an early Neolithic occupation of Malta in in a sediment core extracted from Salina Bay in northeast Malta, while accounting for the high-energy depositional environment and chronological uncertainty associated with radiocarbon dates used to produce associated age–depth models. Each of the analyses was conducted in R and is fully replicable, with scripts, data and outputs contained in a GitHub repository along with further replication instructions (https://github.com/wccarleton/mesoneomalta).

First, we used a standard archaeological phase model to determine start and end boundaries for major depositional phases identified at Latnija. For this model, the excavation team constructed a general Harris matrix relating different contexts to major phases of sediment deposition and artefact accumulation. Thirty-three radiocarbon samples—charcoal from short-lived local shrubs and one bone—recovered from these units were then dated and the dates were placed into an OxCal phase model to estimate phase boundary distributions. All phase boundaries were of the ‘sigma’ type. This boundary allows the tails of the distribution of events (dates) making up abutting phases to overlap. The flexibility reflects the sedimentary fuzziness inherent in the physical boundaries between depositional units. Following previously published guidance⁷², we included a general outlier model along with the phases, allowing for the model to identify potential outliers (events with extreme dates relative to both their phases and the structure of the model as a whole). The modelling identified no significant outliers among the radiocarbon-dated samples given the boundaries we used, as indicated by the posterior probabilities associated with the outlier model that indicate the probability that a given sample is an outlier in the model context (all were 8% or less; most were 4%; Supplementary Table 6).

Next, we used a cleaned regional database of radiocarbon dates associated with securely identified Mesolithic and Neolithic sites or site components from Italy, Sicily, Corsica, Sardinia and Malta. We divided the dates by region and cultural association. Then, we used a simple OxCal phase model to estimate when the Mesolithic phase ended and the Neolithic phase began in each of the regions (details in Supplementary Information 4). We used a single phase for the Mesolithic and one for the Neolithic within each region bookended by flexible ‘sigma’ boundaries, meaning that events at the end of the Mesolithic phase could occur after the estimated start-date boundary for the Neolithic. This flexibility reflects the fact that both phases refer to cultural traditions or packages that are known to have overlapped in space and time throughout the Mediterranean and that have well-established spatio-temporal trends.

Finally, we re-examined the published age–depth model for the Salina Deep sediment core. The core was argued to contain evidence for an early Neolithic in Malta, because it contains findings such as the pollen of domesticated cereals, which was estimated to date to around 8 ka on the basis of an age–depth model. However, the age–depth model used (Bchron), like many sophisticated sedimentation models, assumes monotonicity in the age–depth relationship, which we argue does not apply in the Salina Deep case. Although monotonicity is typically a good working assumption in low-energy depositional environments without evidence of disturbance, Salina Bay in the past and present is a high-energy littoral and fluvial environment that is subject to frequent storms. The core itself contains evidence of marine ingression and many of the radiocarbon dates indicate substantial sediment redeposition, with very old dates near the surface and segments showing a wide radiocarbon temporal spread. Together, this evidence suggests that monotonicity is a poor assumption for Salina Deep and, consequently, that the published age–depth model is overly (unduly) precise because it cannot account for the wide variance in radiocarbon sample dates for many of the core’s segments. To account for this, and produce a model that is more representative of the empirical temporal variance, we used a linear Bayesian regression to model the age–depth relationship. The model recognizes a general relationship in the available age–depth observations indicating a trend toward older dates correlated with depth. However, it also does not assume strict monotonicity, instead focusing on the broad age–depth relationship. We used a custom distribution (based on standard radiocarbon-date calibration) to add a measurement uncertainty component to the model, representing radiocarbon dating and calibration uncertainties. We also used Bayesian imputation to model dates with full posterior uncertainty for a sequence of undated sediment depths (see Supplementary Information 1 for further details).

Zooarchaeological methods

During the 2021 and 2022 field seasons, faunal remains greater than 20 mm in length were piece plotted using a total station, given a unique identifier and bagged. Smaller bone fragments, shells and other faunal remains were recovered through various methods, including an exhaustive programme of wet sieving, flotation and manual inspection of 8-mm, 4-mm, 2-mm, 1-mm and 0.5-mm sieved sinks under microscopy. Here we present a preliminary taxonomic and taphonomic analysis of this faunal material, but note that a full detailed analysis is currently underway that comprises all remains recovered during excavation.

Bones were identified to skeletal element and, for the most part, to broad taxonomic categories (for example, fish and birds), facilitated by relevant literature^73,74,75,76, online resources and comparative material housed at the University of Malta. The taphonomic analysis focused on identifying bone fractures and surface modifications, such as burning, butchery marks (such as cut marks) and carnivore damage (for example, gnawing) following standard protocols^77,78,79,80. Remains are reported as the number of specimens (NRSP) and number of identified specimens (NISP), following a previous report⁸¹. NRSP includes all skeletal remains (bones and teeth) included in this study, whereas NISP is defined as all skeletal elements (bones and teeth) identified minimally to class.

In addition to the piece-plotted bone, we also report here the complete counts of marine fauna for two excavation squares (L2 and N2), reflecting material that was directly recovered and bagged during excavation and material from wet sieving and flotation. Given the very different sediment volumes exposed for the different phases, we chose here to focus at first on these two squares, which offer a good sequence through the phases, to showcase the marine component at the site.

Isotope methods

Nineteen samples, comprising 12 wood mouse (Apodemus sylvaticus) and 7 red deer (Cervus elaphus), were selected for δ¹³C and δ¹⁸O isotope analysis of tooth enamel (Supplementary Table 7). For red deer, molar teeth were targeted for analysis, although the sample set does include one red deer premolar tooth. It should also be noted that because some of these samples are non-overlapping teeth, it is possible that some pseudo-sampling (sampling from the same individual) took place. For wood mouse, whole molar and incisor teeth were used to ensure that the minimum sample size for stable isotope analysis was met.

Before sampling, red deer were cleaned through gentle abrasion with a diamond-tipped drill to remove any adhering material. After cleaning, the same approach was used to sample the tooth enamel along the full length of the buccal surface to ensure a representative measurement for the period of tooth formation. For wood mouse, as much of the dentine was removed as possible using a drill before the remaining whole teeth were crushed using a mortar and pestle, with cleaning of the mortar and pestle using 70% ethanol between samples.

To remove organic or secondary carbonate contaminates, all samples underwent pretreatment, which involved soaking in 0.1 M acetic acid for 10 min followed by three rinses in purified water^82,83. After reaction with 100% phosphoric acid, gasses were analysed using a Thermo GasBench II connected to a Thermo Delta V Advantage mass spectrometer housed at the Department of Archaeology at the Max Planck Institute of Geoanthropology. Carbon and oxygen isotopes are reported as the ratio of heavier to lighter isotopes (¹³C/¹²C or ¹⁸O/¹⁶O) in parts per million (‰) relative to international standards (Vienna Peedee Belemnite, VPDB). δ¹³C and δ¹⁸O values were normalized using a three-point calibration against the international standards IAEA-603 (δ¹³C = 2.5‰, δ¹⁸O = −2.4‰), IAEA-CO-8 (δ¹³C = −5.8‰, δ¹⁸O = −22.7‰) and IAEA NBS 18 (δ¹³C = −5.014‰, δ¹⁸O = 23.2‰), as well as the in-house standard of USGS44 (δ¹³C = −42.2‰).

Reporting summary

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