Here, we report on the vertebrate faunal remains from Phase V deposits recovered during the 2021 to 2024 Latnija excavations. All remains larger than 20 mm in length were piece-plotted using a Leica total station and their coordinates recorded in the Universal Transverse Mercator (UTM) zone 33N (EPSG:32633) geographic coordinate system. To ensure the recovery of even the smallest specimens, all excavated sediments were collected, floated, and the heavy fraction wet-sieved using superimposed eight- to 0.5-mm meshes¹. Taxonomic and anatomical identifications were carried out with the help of modern zooarchaeological comparative material housed at the Department of Classics and Archaeology at the University of Malta, numerous osteological atlases and zooarchaeological guides (e.g., for cervids⁵⁰; for carnivores⁵¹; for tortoises^52,53; for birds^54–56; among others), online resources, and unpublished images. We also used reference collections of Naturalis Biodiversity Center and the zooarchaeology lab of Leiden University (The Netherlands) for the study of avian remains, as well as the National Museum of Natural History of Malta, namely for Anatidae. Skeletal elements and portions were recorded following the coding system of Stiner⁵⁷, except for tortoise’s shell (plastron/carapace) elements, for which we used Thompson and Henshilwood’s nomenclature⁵⁸.

Our main quantification units are the number of identified specimens (NISP) and number of specimens (NSP). NISP includes bone and tooth fragments identified to the lowest taxonomic level possible (e.g., species, genus, family; after⁵⁹). However, here we also utilize NISP values to refer to those remains that were assigned to more general categories (order or class) or body size groups (e.g., small or medium birds); these specimens either lack the diagnostic features that aid specific taxonomic identifications, or belong to very diverse taxonomic groups, such as birds or fish. NSP counts, instead, refer to all the faunal remains reported in this study, including unidentified fragments⁵⁹. In total, we recorded 1,914 NSP and 1,493 NISP. This sample encompasses all the faunal remains that were piece-plotted with the total station during the excavations and recovered as single finds, along with the identified specimens sorted from the heavy fractions recovered by wet-sieving.

A total of 315 specimens were selected for analysis by Zooarchaeology by Mass Spectrometry (ZooMS). The majority (n = 280) derived from the 2 mm screened fraction of the assemblage, which had been recovered through wet sieving followed by flotation. Selection prioritized dense cortical bone fragments, although some spongy bones were included as well. The remaining specimens (n = 35) comprise piece-plotted bones chosen either for their diagnostic potential or because they were of special interest—for example, fish remains—to help refine our understanding of species richness at the site.

All samples were analyzed following an acid insoluble protocol^60–62. For 233 of the samples, approximately 20–50 mg of bone powder was collected using an ultrapure water- and ethanol-cleaned mortar and pestle. Where substantial sediment crust adhered to a specimen, the sampling area was first drilled with an ethanol-cleaned diamond drill bit to remove the sediment. Samples were then demineralized in 500 µL of 0.6 M hydrochloric acid (HCl) at 4°C for 4 hours. After centrifuging, the acid supernatant was removed, and the samples were washed three times with 200 µl of 50 mM Ammonium Bicarbonate (AmBic) solution. The samples were then gelatinized in 100 µl AmBic at 65°C for one hour. Following this, 50 µl of the solution was transferred to a new Eppendorf tube, with 1 µl of trypsin solution (0.4 µg/µl Pierce™ Trypsin Protease, Thermo Scientific), and incubated overnight (~ 16 hours) at 37°C. After incubation, 1 µl of 5% trifluoroacetic acid (TFA) was added to stop tryptic digestion. The samples were purified and desalted using a C18 Ziptip (Piece™ C18 Tips, Thermo Scientific) and were eluted into 10 µl of conditioning solution (0.1% TFA in 50% acetonitrile). A 0.5 µl aliquot was spotted in triplicate onto a MALDI 384 ground steel target plate (Bruker Daltonics) and mixed with 0.5 µl of matrix solution (10 mg/ml α-Cyano-4-hydroxycinnamic acid). For a small subset (n = 82; sample numbers between LAT592–626, LAT674–720), whole bone fragments were used instead of bone powder. For these samples, HCl demineralization was extended to five days, followed by two washes with 200 µl AmBic, two washes with 200 µl of 0.1 M sodium hydroxide (NaOH), and two final washes with 200 µl AmBic.

Samples were analysed at the Institute for Biomedicine and Glycomics, Griffith University, Australia, using a Bruker rapifleX (Bruker Daltonics) in reflector mode, positive polarity, with a mass-to-charge range 800–3500 m/z, laser intensity between 20–60%, and 2000 shots per sample. External mass calibration was performed using a peptide calibration standard (#8206195, Bruker Daltonics) containing a mixture of seven known peptides.

Raw MALDI spectra were first converted to text and mzXML formats in MSConvert (ProteoWizard). Triplicate spectra for each specimen were then processed in R (v4.3.0; ref.⁶³) following⁶⁴ using the ‘MALDIquant’ (v1.22.3) and ‘MALDIquantForeign’ (v0.14.1) packages⁶⁵. This included smoothing with a moving-average algorithm (half-window size = 2), baseline removal using the TapHat method (half-window size = 14), and spectra alignment using the SuperSmoother method (half-window size = 7, signal-to-noise cut off = 3). Spectra were then merged and the baseline corrected once more using the TopHat method. The processed data were exported as .msd files and analysed in mMass⁶⁶, where peaks were re-picked using a signal-to-noise cut-off of 3.5. Taxonomic identifications were made by matching observed peptide masses against a custom reference dataset of Mediterranean taxa^{13,14,60,67,68}.

All faunal remains were examined for bone surface modifications and fracture types to determine the agent(s) of bone accumulation or modification at the cave, as well as post-depositional processes that may have impacted the integrity of the zooarchaeological assemblages. Each specimen was analysed with a 10x hand lens, and when higher resolution or magnifications were required, a binocular microscope. At Latnija, post-depositional modifications include weathering stage, exfoliation, chemical weathering/corrosion, rounding, and trampling (sensu^15,18), as well as the presence of calcitic crusts. Mechanical cleaning of these crusts (after⁶⁹) was attempted but proved to be very time-consuming and often ineffective. As a result, we included in our taphonomic analysis the degree of cortical readability (1: < 50%; 2: 50–75%; and 3: > 75% visible) to establish to which extent bone surface modifications, regardless of the type and accumulating agent, may have been obscured by sediment crusts. Carnivore bite marks, rodent gnawing, and anthropogenic modifications (e.g., burning, percussion damage, peeling, and cut marks) were recorded following^18,70,71, and references therein. Close-up photographs of bone surface modifications were obtained with a Hirox KH8700 digital microscope at the Catalan Institute of Human Paleoecology and Social Evolution (IPHES-CERCA, Tarragona, Spain). Following⁷² and references therein, bone fracture attributes were analyzed to distinguish green fractures—typically linked to human or carnivore activities and marrow extraction—from dry fractures, which tend to result from post-depositional processes like sediment compaction and rock fall. We recorded fracture angles, outlines, edges, and degree of completeness of shaft circumferences for all the mammalian long bone shaft fragments, including proximal and distal ends of each specimen.

Spatial analysis involved examination of the Latnija Phase V faunal dataset, which consisted of 1015 specimens (NSP) piece-plotted during excavations in Trench 4 and included coordinates (XYZ) from the excavation and attribute data from the faunal analysis. The dataset was separated by taxonomic groups (mammals/birds) to assess any spatial patterns in the distribution of fauna related to discard activities, and then by taphonomic indicators (the presence or absence of calcitic crusts covering bone surfaces, and evidence of burning) to investigate patterns in spatial organisation and/or post-depositional alteration. This produced six separate datasets: birds (n = 98), mammals (n = 906), crusts (n = 592), no crusts (n = 294), burnt (n = 302), and unburnt bones (n = 586). During analysis of the Phase V faunal assemblage, a series of refitting sets were also identified (n = 20), and a separate dataset was created to summarise the spatial distribution of these refit sets.

All spatial analyses were performed in R (v4.3.0; ref.⁶³) using the ‘spatstats’ package (v3.0-6; ref.^73,74). Polygons outlining the excavation trench, grid squares, and the cave were constructed using QGIS (v3.30; ref.⁷⁵) and saved as ESRI Shapefiles before being imported into R using the ‘sf’ package (v1.0-14; ref.^76,77). Refits were connected by generating polyline features and horizontal and vertical distances measured between refitting specimens and summarised. A full list of packages and the code used in this study is provided as Supplementary Code.

For analysis of the faunal datasets, spatial data were first checked for duplicates and then converted to point pattern (ppp) objects using the ‘spatstats’ package. Polygons representing the study area and excavation grids were imported and used as spatial windows for the analysis. To assess the spatial distribution of the faunal remains across the Phase V excavation area, kernel density estimation (KDE) was performed using the likelihood cross-validation method to calculate the optimal bandwidth for each dataset. Quadrat counts were calculated to further assess the spatial structure of the faunal assemblage, by subdividing the study area into regular quadrats and counting the number of points within each. To identify significant ‘hotspots’ (i.e., clusters) within each of the datasets, likelihood ratio tests were conducted using the ‘scanLRTS’ function to evaluate and identify locations and spatial scales at which clustering is unlikely to have arisen by chance. To identify areas of clustering, a significance threshold of p < 0.001 was used.

The final stage of the spatial analysis involved testing the point pattern datasets against three different models; the first being that they are distributed randomly (i.e., complete spatial randomness; CSR), the second that there is a positive relationship between them (i.e., they cluster, with points tending to be close together), or thirdly that there is a negative relationship between points (i.e., inhibition, with points tending to avoid each other). The Kinhom and Lscaled functions from the ‘spatstat’ were used to test for inhomogeneous spatial structure. These functions are designed to assess the spatial interaction of points while correcting for varying intensity across the study area. The Kinhom function, or inhomogeneous K-function, is sensitive to clustering or inhibition in point patterns, while the Lscaled function provides a scaled version of the L-function, which normalizes the distance-dependent results. For both functions, envelopes were generated by simulating random point patterns under the null hypothesis of spatial randomness (using an inhomogeneous Poisson process). A total of 99 permutations were performed to generate the envelopes, and the observed function was compared to the simulated envelopes to determine if there were significant deviations from randomness. Maximum Absolute Deviation (MAD) and Diggle–Cressie–Loosmore–Ford (DCLF) tests were then performed on the Kinhom and Lscaled functions to assess whether the observed point patterns were consistent with the null model (Complete Spatial Randomness; CSR), or whether they showed significant deviation from the CSR model due to clustering, inhibition/spatial regularity, or other interaction effects.

Next, possible correlations between the spatial distribution of fauna types (mammals vs aves) and taphonomic indicators (calcitic crusts vs no crusts and burnt vs unburnt) were investigated. This was done using pairwise cross-type K- and L-functions to test for significant interaction between these paired datasets, with envelopes again generated using 99 permutations. Finally, quadrat count chi-square tests were performed at different resolutions (5x5 m, 8x8 m, 10x10 m, 12x12 m) to test whether the spatial distribution of each type deviated from complete spatial randomness (CSR) at different scales, producing a chi-squared statistic (X²), along with the corresponding degrees of freedom. For all statistical tests the results were assessed based on p-values, with evidence considered as very strong (p ≤ 0.001), strong (p ≤ 0.01), moderate (p ≤ 0.05), or non-significant (p > 0.05) following the approach advocated by⁷⁸.

In order to evaluate ungulate skeletal element representation and constrain simultaneously the character of initial accumulation related to a given transport strategy/depositional process, and the subsequent occurrence of bone density-mediated attrition, we applied a Bayesian method proposed by¹⁶. First, minimum number of elements (MNE) (sensu⁵⁹), were derived from NISP by considering all the complete and fragmentary specimens and quantifying the minimum number of skeletal portions or unique bone landmarks for a given element. Second, Minimum Animal Units (MAU) were estimated following^79,80 by dividing the observed MNE values by the expected MNE for each element in a complete cervid skeleton. Then, we used %MAU and the ‘BaskePro’ package in R (v.1.1.1; ref.¹⁷). Based on a Monte Carlo Markov Chain sampling, this Bayesian method considers two parameters: 1) alpha (α), which informs about transport/skeletal element selection, and its value can range from − 1 (the original accumulation was dominated by axial elements) to 1 (mostly appendicular elements); and 2) beta (β) or degree of attrition, which correlates the survivorship of bone elements to their maximum mineral density, and it can take any value between 0 (no attrition) to 10 (maximum attrition)^16,17.

Additionally, log-normal family generalized linear models were fitted to our data to model: 1) the effects of bone mineral density on percent survivorship^81,82; and: 2) %MAU as a function of various utility indices, including the Modified General Utility Index (MGUI)⁸³, Simplified Meat Utility Index (SMUI)⁸⁴, Food Utility Index (FUI)⁸⁴, Corrected Food Utility Index (CFUI)⁸⁵, and Unsaturated Marrow Index (UMI)⁸⁶. We utilized bone density values from American deer (Odocoileus spp.)^81,82, and utility values calculated for caribou/reindeer (Rangifer tarandus)^83–86 since they are the closest taxa to Mediterranean red deer (Cervus elaphus) for which these data are available. Following⁸⁷, statistical significance and goodness of fit for all log-normal family generalized linear models were examined with a likelihood ratio test in the ‘lmtest’ package (v. 0.9–40; ref.⁸⁸) and D² statistic (the proportion of variation explained by the models) in the ‘modEva’ package (v.3.40; ref.⁸⁹).

Considering the abundant avian remains at the site, skeletal element profiles to distinguish anthropogenic and non-anthropogenic bird assemblages were explored, focusing on the wing-to-leg ratio (after⁹⁰; and modified by⁹¹) and proximal-to-distal ratio (following⁹²). It has been shown that these ratios differ between bird assemblages accumulated by different predators or resulting from in situ mortality (e.g.^21,90,93), where they can also be used to assess the influence of bone density-mediated attrition in structuring element abundances (e.g.^94,95; but see also⁹¹). In this study, the wing-to-leg and proximal-to-distal ratios are based on MNE values, which were estimated considering all avian remains assigned to the medium body size group, since they are the most numerous at the site. To have a more robust sample, we also calculated them for all the different bird body size categories combined.

Data on dental eruption and wear stages for red deer specimens (after²⁷, for the mandibular teeth; and after⁹⁶; for the upper teeth), as well as long bone epiphyseal fusion for cervid/ungulate remains were collected to reconstruct mortality profiles. Following²⁷, the different eruption and wear stages were grouped into three age categories: 1) juveniles (only deciduous teeth), 2) prime-aged adults (permanent teeth with no occlusal wear up to medium wear), and 3) old adults (permanent teeth with medium-advanced to advanced wear). To avoid pseudo-replication or counting the same individual more than one, aging of dental remains should focus exclusively on lower deciduous premolars (Dp₄) and lower permanent fourth premolars (P₄; ref.²⁷) or third molars (M₃; ref.⁹⁷). When combined, these sets of elements represent the complete lifetime of an individual since the deciduous tooth is shed and replaced by its permanent counterparts²⁷. Yet, due to the available sample size, as a second dataset, all the dental lower and upper premolars and molars were also combined to provide a more robust sample. We plotted the proportion of juveniles, prime, and old cervids on a triangular diagram with 95% confidence interval using Weaver et al.’s cross-platform computer program⁹⁸. Our results were then superimposed on the diagram proposed by⁹⁹ for red deer, whose different zones have been mathematically modeled considering the age patterns of living populations.

Metric analyses of deer and ungulate specimens follow¹⁰⁰, and all linear measurements were taken using electronic digital callipers. Due to the high degree of fragmentation affecting the bones from Latnija and limitations associated with small sample sizes, the Size Variation Index (SVI) was used to examine body size differences and enable comparisons between early Holocene Maltese deer populations and their Late Pleistocene and Holocene Sicilian and southern Italian counterparts. Like other scaling techniques (e.g.^101,102), the SVI standardizes linear measurements of different skeletal elements based on a reference sample, which allows us to provide a more robust sample and make comparisons between assemblages for which different measurements or sets of elements are available¹⁰³. However, only one type of measurement per specimen was included in our analysis to avoid pseudo-sampling (Supplementary Table 22). Our standard sample (mean [M] and standard deviations [SD] values) was calculated using metric data provided by¹⁰⁴ for Late Pleistocene deer remains recovered in Sicily due to its geographical proximity; yet, we only considered those measurements that were consistent with von der Driesch’s standards¹⁰⁰ (Supplementary Table 23). The SVI values were calculated as (X - M)/SD, where X is a given measurement of a bone, and M and SD are the mean and standard deviation values in the standard sample (after¹⁰³). Although Di Stefano et al.¹⁰³ listed several southern Italian deer assemblages with metric data, we only considered those for which single values (rather than averages and standard deviations of a set of measurements) were available in the literature (i.e.^105–110). Metric data for the Sicilian deer are from¹⁰⁴. To compare the median SVI values of the Maltese, Sicilian, and southern Italian deer, we performed a Kruskal-Wallis rank sum (non-parametric) test in R. Since the results showed strong evidence that median SVI values are different between groups, the analysis of body size was followed by a pairwise comparison to identify which deer assemblages differ between each other.

Because prey body mass not only allows us to make inferences about past capacity and meat availability, but also affects interpretations of skeletal element representation in relation to hunting grounds, transport patterns and body part selection by hunter-gatherers (e.g.¹¹¹), an allometric based method for estimating body mass was applied to selected deer postcranial specimens from Latnija. The predictive equation is as follows: Log body mass (Kg) = b(log X) + a, where b is the slope, X is one of our linear measurements (in mm), and a is the intercept of the regression (sensu^112,113). a and b values for cervids are provided by¹¹³. Our estimations focus exclusively on the diameters of proximal and distal epiphyses (transverse sections) of long bones, which tend to be more highly correlated with body mass than bone length measurements (e.g.¹¹³). The elements used for body mass estimations in our study are humerus (H4, H5), radius (R4, R5), metacarpal (Mc2, Mc4), tibia (T4), and metatarsal (Mt2, Mt4) (ref.¹¹²), which are also equivalent to von den Driesch’s measurements¹⁰⁰ (Table 2 and Supplementary Table 22).