An infant burial from Arma Veirana in northwestern Italy provides insights into funerary practices and female personhood in early Mesolithic Europe

Fieldwork at Arma Veirana and all sampling of the AVH-1 skeleton and other artifacts was conducted under the auspices of Soprintendenza Archeologia, Belle Arti e Paesaggio per la città metropolitana di Genova e le province di Imperia, La Spezia e Savona via a permit formally issued to coauthor Fabio Negrino, following all relevant requirements of the permitting authority.

Fieldwork and site documentation

Spatial relationships in the Arma Veirana excavation have been documented using total stations and photogrammetry to establish 3D coordinates for all artifacts recovered in situ (no size cutoff). Details concerning establishment of the excavation grid using control points (Fig. S1; Table S2) and other details are provided in the Supplementary Information. Standard field-based descriptions of deposits and stratigraphy were combined with micromorphology. Micromorphological analysis of intact blocks of sediment allowed for high resolution documentation of the burial’s geological context, including the composition and origin of the burial fill (see Supplementary Methods).

Excavation of the delicate AVH-1 burial was not conducted en bloc; due to the surroundings’ rugged terrain, transporting a burial block to the lab was not possible. Instead, each piece was removed on site. To document the 3D spatial relationships of skeletal remains and grave goods to one another, total station data were combined with a system of progressive photogrammetry (Fig S2). Further technical details are provided in the Supplementary Methods.

Radiocarbon dating and carbon stable isotopes

With permission of the Soprintendenza Archeologia, Belle Arti e Paesaggio per la città metropolitana di Genova e le province di Imperia, La Spezia e Savona, and following all relevant guidelines issued by that authority, one human vertebral neural arch from the Arma Veirana burial (plotted find #9237; Fig. S3) was pretreated at the Department of Human Evolution at the Max Planck Institute for Evolutionary Anthropology, Leipzig, Germany, using the method described both in39,40 for bone samples 80 µm) particles. The gelatine was ultrafiltered41 with Sartorius ‘VivaspinTurbo’ 30 KDa ultrafilters. Prior to use, the filter is cleaned to remove carbon containing humectants42. The samples were lyophilized for 48 h.

The collagen was weighed into a pre-cleaned tin cup and sent to the Curt-Engelhorn-Centre for Archaeometry Klaus-Tschira-AMS facility in Mannheim, Germany (lab code: MAMS) for graphitization and dating with the MICADAS-AMS43. To monitor contamination introduced during the pre-treatment stage, a sample from a cave bear bone, kindly provided by D. Döppes (MAMS, Germany), was extracted along with the batch AV144. To assess the preservation of the collagen yield, C:N ratios, together with isotopic values must be evaluated. The C:N ratio should be between 2.9 and 3.6 and the collagen yield not less than 1% of the weight45. Stable isotopic analysis was conducted at MPI-EVA, Leipzig (Lab Code R-EVA), using a ThermoFinnigan Flash EA coupled to a Delta V isotope ratio mass spectrometer.

Several other samples collected in 2017 at the beginning of the burial excavation were dated at the University of Oxford Radiocarbon Acceleration Unit (ORAU) to provide broader contextual information to inform interpretations of the burial and site formation. These included faunal bone and charcoal samples that derive from the burial fill or adjacent sediments. The ORAU protocols include the correction of isotopic fractionation using δ13C values measured using the accelerator mass spectrometer, which are also reported independently using results from a stable isotope mass spectrometer (to ± 0.3 per mil relative to VPDB). Details of the ORAU protocols can be found elsewhere46,47.

Dental analyses

Virtual histology

Human teeth permanently record growth information during their formation allowing assessment of tooth formation times, stresses experienced during development, and age-at-death (in infants). Dental ontogenetic studies rely on the rhythmical growth of enamel and dentine, producing short and long period incremental markings, visible in longitudinal thin sections of teeth48. The deposition of enamel and dentine matrix is subject to inner biological rhythms: a circadian growth process that produces prismatic cross striations (CS) in enamel and von Ebner’s lines in dentine, and a longer periodicity represented by Retzius lines in enamel and Andresen lines in dentine16,49. Stress events strong enough to disrupt development leave marks in the corresponding position of the enamel or dentine front, visible as Accentuated Lines (ALs)49. Usually, the first accentuated line, characterizing all the deciduous teeth and the first permanent molar of individuals that survive the perinatal stage, is the Neonatal Line (NL), which separates the tissues formed prenatally from that growing after birth50. ALs in the prenatal enamel are relatively rare and indicate physiological stresses affecting the mother or fetus51. The NL is used to calibrate the chronology of stresses and the time of pre- and postnatal crown growth. Histomorphometry of dental enamel allows the collection of parameters such as the Daily Secretion Rate (DSR, i.e. the speed at which the ameloblast—the enamel forming cells—moves towards the outer surface of the tooth) and, for still growing crowns, the age-at-death as the time spent from the NL to the end of enamel formation16,51,52.

Virtual histomorphometry of three deciduous tooth crowns (AVH-1d, left upper central deciduous incisor; AVH-1c, left upper lateral deciduous incisor; AVH-1g, right upper deciduous first molar) was conducted to identify the NL, assess chronological age-at-death, and investigate biological life history. Precautions were adopted to preserve ancient DNA53. The dental deciduous crowns were imaged through synchrotron radiation computed microtomography (SRμCT)17,51,54 at the SYRMEP beamline of the Elettra Sincrotrone Trieste laboratory in Basovizza (Trieste, Italy)55. Details of the scanning protocol are provided in the supplementary materials.

Virtual histological sections passing through the bucco-lingual plane at the tip of the dentine horn of the incisors and of the mesio-lingual cusp (protocone) of AVH-1g were derived from the synchrotron radiation computed microtomography volumes (Figs. 2A, S4, S5). Each virtual slice was obtained using the average Z projection of 6 to 10 consecutive reformatted slices using the Z-Project tool of the software Fiji, thus resulting into 9 to 30 μm thick reformatted slices (Figs. S4 and S5).

Proteomic analysis: amelogenin-based sex estimation

Proteomic analyses through LC–MS/MS were conducted on the enamel of the left maxillary first molar to estimate the sex of AVH-1. Before sampling, laboratory-based X-ray microCT scan of the tooth was undertaken by use of a desktop Skyscan 1072 system (Bruker Corp., Kontich, Belgium). Scanning parameters were 50 kV and 197 µA for Voltage and current of the X-ray source, 1 mm aluminum filter, an exposure time of 5936 ms, image averaged on 2 frames, a total scan angle of 180° with an angular step of 0.9°. Reconstruction was obtained with NRecon software (Bruker Corp., Kontich, Belgium), with an isotropic voxel size of 10.42 µm, beam hardening correction (%) = 25, and ring artifact correction = 1.

Amelogenin-based sex estimation56enables the confident identification of the sex of an individual even when DNA seems compromised22. The analytical protocol is described in Lugli et al.23. A small fragment of the tooth (~ 10 mg; including both dentine and enamel) was sonicated with MilliQ water and quickly pre-cleaned with 5% HCl. Then, the specimen was digested through 250 μl of 5% HCl (suprapur grade). This solution was desalted and purified using a C 18 silica SpinTip (Thermo Scientific) and dried down overnight under a class 100 laminar-flow hood in the clean room facility of the Department of Chemical and Geological Sciences (University of Modena and Reggio Emilia). Dry peptides were re-suspended using 35 μl of a water:acetonitrile:formic acid mixture (95:3:2) and analysed by LC–MS/MS (Dionex Ultimate 3000 UHPLC coupled to a high-resolution Q Exactive mass spectrometer; Thermo Scientific). A run time of 90 min was employed for the specimen and the blanks; see23 for details. The analysis of the AVH-1’s tooth was duplicated. To check the presence of amelogenin proteins we searched the ion chromatograms22,23, focusing on specific peptides as SM(ox)IRPPY (AMELY; [M + 2H]+2 440.2233 m/z) and SIRPPYPSY (AMELX; [M + 2H]+2 540.2796 m/z). Additionally, raw data were converted in Mascot generic format and searched against UniProt (constrained to Homo sapiens) and cRAP (contaminant database, 116 sequences). No proteolytic enzyme was selected in search parameters and one missed cleavage allowed. Deamidated asparagines/glutamine (NQ) and oxidated methionine (M) were set as variable modifications. Error tolerance was set at 10 ppm for the precursor ions and 0.05 Da for the product ions. False discovery rate was estimated through an automatic decoy database, with a probability threshold trimmed to a FDR < 1%. A specific protein was considered as identified if at least two significant peptides were observed.

Ancient DNA extraction, library preparation, mitochondrial DNA capture and sequencing

After removing a thin layer of surface from the vertebral fragment of AVH-1 (plotted find #9237; Fig. S3), we drilled 9.5 mg and 8.1 mg of bone powder with a sterile dentistry drill. We produced 500 µl of lysate from each sample and used 150 µl aliquots for automated silica-based DNA extraction57, choosing binding buffer option ‘D’. Both DNA extracts were converted into single-stranded DNA libraries, which were subsequently quantified, amplified, and labelled with two unique indices58. Negative controls for the DNA extraction and the library preparation were carried along through all steps. The amplified libraries were pooled and sequenced on an Illumina HiSeq 2500 in paired end mode (2 × 75 cycles) with two index reads59. An aliquot of each library was further enriched for human mitochondrial DNA (mtDNA) in two successive rounds of hybridization capture60,61,62. Enriched libraries were pooled and sequenced on an Illumina MiSeq.

Paired-end sequence reads were overlap-merged into single-molecule sequences and adapters were trimmed using leeHom63. Merged sequences were mapped using the Burrows-Wheeler Aligner (BWA)64 with ancient DNA parameters (‘–n 0.01 –o 2 –l 16500’)65, either to the modified human reference GRCh37 from the 1000 Genomes project (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/) for the shotgun data, or to the revised Cambridge Reference Sequence (rCRS) for the mtDNA capture data. All downstream analyses were restricted to sequences with indices that perfectly matched the expected index combinations. We used bam-rmdup (version: 0.6.3; https://github.com/mpieva/biohazard-tools) to remove PCR duplicates and SAMtools (version: 1.3.1)66 to filter for fragments that were longer than 35 base pairs (bp) and had a mapping quality greater or equal to 25. Tables S2 and S4 summarize the number of DNA fragments retained after each step.

To evaluate if some of these fragments stem from authentic ancient DNA, we determined the frequency at which cytosines (C) are substituted by thymines (T) at their ends67. The elevated frequencies of C-to-T substitutions in the libraries of AVH-1 (Table S2 and S5) indicate the preservation of ancient DNA molecules in the specimen. Furthermore, these frequencies remained stable after filtering for fragments with a C-to-T substitution at the opposing end (‘conditional’ substitutions) (Tables S3 and S5), indicating that the majority of the data come from authentic ancient DNA68.

Reconstruction of the mitochondrial genome and phylogenetic analysis

Using schmutzi 69 (parameters: ‘—notusepredC —uselength’), we estimated 1% (95% confidence intervals (CI): 0–2.0%) of present-day human DNA contamination among all mtDNA fragments obtained from the AVH-1 libraries (Table S4). We next reconstructed the full mtDNA genome of AVH-1, first using all mapped fragments longer than 35 base pairs with a mapping quality of at least 25, and second using only fragments with a C-to-T difference to the reference genome at the first three and/or last three terminal positions. We called a consensus base at each position along the mtDNA genome that was covered by at least three DNA fragments and where at least 67% of fragments carried an identical base, as detailed in68. The reconstructed mitochondrial genomes from all fragments and from putatively deaminated fragments were identical, and identical to the consensus sequence called using schmutzi. We identified a haplogroup of AVH-1 to be U5b2b using HaploGrep270,71 and PhyloTree database (PhyloTree.org, build 17)72.

A tip date for the AVH-1 mtDNA was estimated using the Bayesian method implemented in Beast2 (version 2.4.8)73. The reconstructed mitochondrial genome of AVH-1 was aligned to the mtDNA genomes of 54 present-day74 and 52 ancient modern humans of known radiocarbon age62,75,76,77,78,79,80,81,82, which served as calibration points for tip dating, using MAFFT v7.27183 (Table S6). The Vindija 33.16 Neandertal mtDNA genome74 was used as an outgroup. Using jModelTest284 we determined the best-fitting substitution model for this dataset to be Tamura-Nei 93 with a fixed fraction of invariable sites and gamma distributed rates (TN93 + I + G). Models of rate variation and tree priors were investigated following79,82. We determined the strict clock model and Bayesian skyline to be the best fit to the data following a marginal likelihood estimation (MLE)30 analysis for model comparison and best support assessment. The resulting tree was visualized using FigTree (version: v1.4.2) (http://tree.bio.ed.ac.uk/software/figtree/).

Sex determination from nuclear DNA

Sex of AVH-1 was determined from the shallow shotgun sequencing by counting the number of fragments which aligned to the X chromosome and the autosomes. The analysis was first performed using all mapped fragments longer than 35 base pairs with a mapping quality of at least 25, and then using only fragments with a C-to-T difference to the reference genome at the first three and/or last three terminal positions. Sex was determined based on the expected ratios of X to (X + autosomal) fragments for male and female individuals85.

Analysis of the Burial Artifacts

Shell beads and pendants

Perforated shells were observed under ~ 20×–150× magnification using a DinoLite. Microscope images were used to take general measurements of the shells and their perforation, as well as to document the location of use-wear and ochre traces. A selection of 60 Columbella rustica, the single Turritella sp., and the four pendants were examined using a stereoscopic microscope (ZEISS AXIO Zoom V16) at the Diet and Ancient Technology (DANTE) Lab of La Sapienza,
https://www.nature.com/articles/s41598-021-02804-z