Whole-genome ancestry of an Old Kingdom Egyptian

Provenance and ethics The human remains were excavated from the Nuwayrat necropolis near Beni Hasan, Egypt. They were donated between 1902 and 1904 by the Egyptian Antiquities Service to the members of the Beni Hasan excavation committee and subsequently donated to the Institute of Archaeology, University of Liverpool and exported under the John Garstang export permit. The human remains were then donated to the World Museum (previously the Liverpool City Museum) in 1950. Sampling permit was granted by the World Museum. Ancient DNA extraction, library preparation and sequencing Sampling and DNA extraction of seven permanent teeth belonging to an individual from Nuwayrat were carried out in dedicated ancient DNA facilities at Liverpool John Moores University. Library preparation and sequencing were carried out at The Francis Crick Institute (Supplementary Table 1). Before subsampling, the teeth were decontaminated by wiping with 1% sodium hypochlorite, followed by wiping with molecular biology grade water and ethanol. Approximately 44–66 mg of cementum-enriched powder was extracted from each tooth using a Dremel drill at the lowest possible rotations per minute (5,000 rpm). DNA was extracted using 1 ml of extraction buffer consisting of 0.45 ml of 0.5 M EDTA (pH 8.0) and 10 μl of 10 mg ml−1 proteinase K per 50 mg of bone powder. The mixture was incubated overnight (approximately 18 h) at 37 °C and purified on the High Pure Viral Nucleic Acid Large Volume Kit (Roche) using a binding buffer described in ref. 63 and QIAGEN buffer PE. DNA was eluted in approximately 100 μl of QIAGEN elution buffer. Extracts were turned into single-stranded DNA libraries39 (without treatment to remove uracils), double-indexed64 and then underwent paired-end sequencing on a HiSeq 4000 to approximately seven million reads per library for initial screening (Supplementary Table 1). All samples were processed alongside negative lysate and extraction controls and positive and negative library controls. On the basis of the assessment of the initial sequencing results, two libraries were selected for extra rounds of deeper sequencing on the NovaSeq 6000 and NovaSeq X platforms, following the selection of fragments greater than 35 bp using polyacrylamide gel electrophoresis65 for the library built from the NUE001b5e1 extract (Supplementary Table 1), with a resulting total of 8.3 billion 2 × 100 sequence pairs. Radiocarbon dating New radiocarbon dating was generated for the individual that yielded DNA from Nuwayrat (NUE001) by Beta Analytic using accelerator mass spectrometry. We directly dated the upper-left third molar (NUE001b3) and lower-left first premolar (NUE001b5), both of which yielded DNA that was deep sequenced (Supplementary Table 1). The results are reported in Supplementary Table 2. The femur of this individual was previously radiocarbon dated66,67 (Supplementary Information section 1 and Supplementary Table 2). All dates were calibrated using OxCal v.4.4.4 (ref. 68) with atmospheric data in IntCal20 (ref. 69). We also combined the three independent dates using the R_Combine() function in OxCal70 (Supplementary Table 2). We rounded the calibrated dates outwards to 10 years unless error terms were smaller than ±25 bp, in which case we rounded outwards to 5 years71. Isotope analysis Dental collagen and enamel were extracted from the lower-left second molar. Dentine collagen was extracted for carbon (δ13C) and nitrogen (δ15N) isotope analysis following a modified Longin method72,73. Mass spectrometry was performed using a Flash 1112 series elemental analyser coupled with a Finnigan DELTA V Advantage (Thermo Fisher Scientific) using established protocols74. Analytical precision (1σ) of the in-house calibrated standards74 were 0.08 and 0.07 for δ13C and δ15N, respectively. For enamel, after surface abrasion, a slice (3.6 mm wide) was extracted, and all adhering dentine was removed. Two fragments were powdered, one of which was pre-ultrasonicated. A minimum of 3.0 mg was analysed for the oxygen isotope composition of enamel carbonate (δ18O c ). Samples were acidified for 5 min with more than 100% ortho-phosphoric acid (density approximately 1.9 g cm–3) at 70 °C and analysed in duplicate using a MAT 253 dual-inlet mass spectrometer (Thermo Fisher Scientific) coupled to a Kiel IV carbonate preparation device using established protocols74. Isotope values are reported as per mille (18O/16O) normalized to the Vienna Pee Dee Belemnite (VPDB) scale using an in-house carbonate standard (BCT63) calibrated against NBS19. The long-term reproducibility for δ18O BCT63 is ±0.04‰ and ±0.03‰ for δ13C (1σ). The oxygen carbonate values (δ18O C VPDB ) were converted to the Vienna Standard Mean Ocean Water (VSMOW) scale75 and phosphate (δ18O P VSMOW )76. The remaining enamel fragment (56.3 mg) was cleaned in an ultrasonic bath, digested in 8 M HNO 3 and heated overnight at 120 °C. Sr-Spec was used for strontium extraction, following the revised version of Font et al.77. Once column-loaded in 1 ml of 8 M HNO 3 , matrix elements were eluted in washes of 8 M HNO 3 , and samples were placed on a hotplate (120 °C) overnight, with a repeat pass following. The sample was redissolved in 2% HNO 3 , and the 87Sr/86S ratio was measured using a Neoma multi-collector inductively coupled plasma mass spectrometry with tandem mass spectrometry (MC-ICP–MS/MS, Thermo Fisher Scientific). Instrumental mass bias was corrected for using the exponential law and a normalization ratio of 8.375209 for 88Sr/86Sr (ref. 78). Residual krypton (Kr) and rubidium (87Rb) interferences were monitored and corrected using 84Kr and 86Kr (83Kr/84Kr = 0.20175 and 83Kr/86Kr = 0.66474; without normalization) and 85Rb (85Rb/87Rb = 2.5926), respectively. The accuracy of the method was assessed by measuring the EC-5 coral standard (87Sr/86Sr: 0.709171 ± 0.000016 (2σ; n = 14), consistent with the expected value for seawater). The data were also corrected against a National Institute of Standards and Technology Standard Reference Material 987 value of 0.710248 (ref. 79). The procedural blank was less than 75 pg of Sr, negligible relative to sample Sr. Osteological analyses Following element inventory, our determination of the Nuwayrat individual’s sex was based on standard morphological indicators across the skeleton (protocol in Buikstra and Ubelaker25). Ageing was estimated from the dentition, cranium and postcrania25,27,80,81,82,83,84. For stature, several approaches were used85,86,87, with the most likely estimate based on direct stature reconstruction of ancient Egyptians following ref. 26. Biological affinity was assessed from two long-recognized methods: dental non-metric traits88 and craniometrics (for example, Howells89). First, the rASUDAS application was accessed (https://osteomics.com/rASUDAS2/)90. It used up to 32 crown and root traits for comparison with data from seven global population samples. Second, the craniometric approach used the CRANID program CR6bIND, with 29 measurements for comparison with a database of 74 premodern through recent global samples, including Late Dynastic Egyptians and ancient West Asians91. Our recording and description of skeletal pathology, related primarily to age-related breakdown, follow accepted methods25,92,93. This and activity-induced musculoskeletal stress markers (details in previous studies28,29,94) were used to ascertain the level of physical activity. Although not without criticism95,96, they have been used to infer occupation by identifying common positions and movements in life. For that purpose, the latter were compared with illustrations of individuals engaged in a range of common jobs, as depicted on ancient Egyptian tomb walls and in statuary (Supplementary Information section 2). Facial reconstruction and depiction Craniofacial analysis and facial reconstruction from skeletal remains were carried out using three-dimensional laser scan data of the skull (collected using an Artec Space Spider scanner), Touch X haptic device and Geomagic Freeform software97. Egyptian male data98 were used to estimate facial tissues at anatomical points across the skull surface. The muscles of the head and neck were imported from the Face Lab database and remodelled to fit the skull following anatomical guidelines99. Morphometric standards were used99,100 to estimate facial feature morphology, such as eye and nasal shape, lip and ear pattern and structural creases. A final facial depiction was produced using two-dimensional photo-editing software. It is important not to consider a facial depiction as a portrait or definitive image because it can only visualize the available information101. In this case, although DNA analysis indicated the most probable population of origin, there was no evidence in relation to skin colour and hair colour. Therefore, the facial depiction was produced in black and white without head hair or facial hair (Supplementary Information section 3). Bioinformatics data processing and authentication Read alignment was performed following the pipeline in the study of Swali et al.102. Samples were processed through the nf-core/eager v.2.3.3 pipeline103. First, adaptors were removed, paired-end reads were merged and bases with a quality below 20 were trimmed using AdapterRemoval v.2.3.1 (ref. 104) with –trimns –trimqualities –collapse –minadapteroverlap 1 and –preserve5p. Merged reads with a minimum length of 35 bp were mapped to the hs37d5 human reference genome with Burrows-Wheeler Aligner (BWA-0.7.17 aln)105 using -l 16500 -n 0.01 -o 2 -t 1 (ref. 106). Duplicate reads were removed using DeDup v.0.12.8 (ref. 107). Finally, we removed the alignments with mapping quality below 30 and containing indels. We used mapDamage v.2 (ref. 108) to visualize the substitution distribution along the reads and evidence the presence of deaminated molecules typical of ancient DNA. Contamination was estimated using three different data sources: (1) genome-wide present-day contamination using the conditional substitution rate109 computed using PMDtools v.0.60 (ref. 110); (2) present-day mitochondrial DNA-based contamination using schmutzi (commit be61017)111; and (3) chromosome X contamination on libraries assigned as male using ANGSD v.0.933 (ref. 112), restricted to the non-recombining region of chromosome X. All the libraries from NUE001 show little to no contamination, except two libraries with sequencing identification numbers SKO719A1706 and SKO719A1709 (Extended Data Fig. 3 and Supplementary Table 1). Molecular sexing The biological sex of the sequenced individual was determined using the R y parameter113, which is the ratio of the number of alignments to the Y chromosome (n y ) to the total number of alignments to both sex chromosomes (n x + n y ), R y = n y /(n x + n y ). All libraries are consistent with NUE001 being karyotypically male, except the results from the library SKO719A1706 consistent with being female, which is probably a result of contamination (Supplementary Table 1). SNP calling in the Nuwayrat individual We merged the sequencing data from five libraries from the Nuwayrat individual showing an absence of present-day human DNA contamination, yielding a total of 135,606,409 mapped unique reads of 44.63 bp on average, resulting in an average genome-wide coverage of 2.02×. We called pseudo-haploid positions using SAMTools v.1.9 mpileup114 with options -B -R -Q30 and SequenceTools 1.5 (ref. 115) with options –randomHaploid and –singleStrandMode. This approach leverages the single-stranded library preparation to computationally remove the effects of cytosine-deamination-derived sequence errors. Specifically, at C/T SNPs, it removes all bases that are aligned onto the forward strand; at G/A SNPs, it removes all bases on that aligned to the reverse strand. This allows for a confident pseudo-haploid genotyping even also at CpG context transitions, which are mostly not repaired by the uracil-DNA glycosylase (UDG) treatment owing to methylation116. Uniparental marker determination We obtained the mitochondrial DNA consensus of the Nuwayrat individual from endogenous reads, removing the bases with quality below 20 (-q 20) using schmutzi111. The mitochondrial haplogroup was assigned using Haplogrep 3 (ref. 117). The chromosome Y haplogroup was obtained using pathPhynder118 with the parameter -m ‘no-filter’, on the basis of approximately 120,000 SNPs extracted from worldwide present-day and ancient male chromosome Y variation and the International Society of Genetic Genealogy v.15.73 (http://www.isogg.org). Comparison dataset We merged the genome of the Nuwayrat individual with a comparison dataset of 977 ancient individuals2,3,4,5,20,43,45,46,47,53,54,55,60,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160 and 4,040 modern individuals43,46,122,141,158,161,162,163,164,165,166,167,168,169 genotyped on either the Human Origins array169 (‘Human Origins’ dataset) or the 1.2 million SNP array (‘1240k’ dataset)139 (Supplementary Table 3). Most genotypes were directly accessed from the Allen Ancient DNA Resource v.54.1 (ref. 170). We added nine ancient genomes from Morocco45 and 13 ancient genomes from Mesopotamia120 from raw mapped Binary Alignment Map (BAM) files processed following the above-mentioned bioinformatic pipeline, with two modifications: (1) for the double-stranded UDG-treated genomes from ref. 45, we trimmed the first and last three bases of the reads and then called pseudo-haploid genotypes at both transition and transversion sites; and (2) for the non-UDG-treated genomes from ref. 120, we called pseudo-haploid genotypes at transversion sites only. We included 100 present-day Egyptian genomes from ref. 164 in both datasets. Individuals related up to the second degree, as detected in previous studies, were excluded. Principal component analysis We computed two PCA on present-day individuals from the ‘Human Origins’ dataset using 593,124 substitutions through SMARTPCA (eigensoft v.6.1.4)169. For the first analysis, we kept 3,233 individuals from across the world and projected NUE001 on the resulting components. For the second PCA, we kept 722 present-day individuals from North Africa, West Asia and the Caucasus and projected NUE001 together with 781 ancient genomes from North Africa, West Asia and the Caucasus. Both analyses used transversions only (111,208 SNPs). ADMIXTURE clustering We used a model-based clustering approach from the program ADMIXTURE v.1.2 (ref. 40) to estimate the ancestry components from genomes in the ‘Human Origins’ dataset. All genomes were transformed into pseudo-haploid sequences, and transitions were removed. The remaining 111,208 positions were subsequently pruned for SNPs in strong linkage disequilibrium using PLINK v.1.9 (ref. 171), with the parameter –indep-pairwise 200 25 0.4 to yield a final set of 71,202 transversion SNPs. ADMIXTURE was run with cross-validation enabled using --cv flag for all ancestral population numbers from K = 3 to K = 20. Runs of homozygosity The presence and length of runs of homozygosity greater than 4 cM in the Nuwayrat genome were estimated using hapROH v.0.64 (ref. 41) on the 1.2 million SNP set of sites. qpAdm modelling For all qpAdm modelling in this study, we estimated the ancestry proportions as a mixture of a set of left (source) rotating populations differentially related to a set of right (outgroup) populations using ADMIXTOOLS 2 (ref. 42) qpadm_rotating() with the option maxmiss = 0.1, removing genotypes missing in more than 10% of populations. We restricted the analysis to genomes with both transitions and transversions (half/plus UDG-treated libraries or single-stranded libraries called using SequenceTools115 --singleStrandMode), removing CpG sites, to increase the robustness of the models. We considered only models with three or less sources. We restricted the fixed set of outgroup to populations distantly related to any left populations and with genomes greater than or equal to 2×: Ju_hoan_North.DG, Ethiopia_4500BP.SG, Latvia_HG_UDG, USA_Ancient_Beringian.SG, Vanuatu_400BP_UDG, Japan_HG_Jomon_UDG and China_NEastAsia_Coastal_EN_UDG. This analysis was conducted on the 1240k dataset. These parameters are always true unless otherwise stated. We ranked the non-rejected models first on the basis of the minimal number of source populations, assuming that a fewer number of source populations is more parsimonious. Then, if several models with the same number of sources are not rejected, we considered the P value, given that the number of SNPs in the rotating models are nearly equal (10% missingness allowed between populations). Supplementary Information section 4 details all models tested. Nuwayrat genome ancestry modelling We first estimated the Nuwayrat genome (NUE001) and contemporary North African and West Asian populations (Levant_BA (also for each of the eight archaeological sites separately), Anatolia_BA and Morocco_MN) ancestry proportions as a combination of distal Neolithic populations from North Africa and West Asia (Morocco_Epipaleolithic, Anatolia_Neolithic, Levant_Neolithic, Zagros_Neolithic and Caucasus_Neolithic). This analysis was carried out on 433,280–558,848 SNPs. Then, we estimated NUE001, Bronze Age Levant (also for each of the eight archaeological sites separately) and Bronze Age Anatolia ancestry components, adding more proximal Neolithic and Chalcolithic North African and West Asian populations as potential sources (Morocco_EN_ktg, Morocco_MN, Anatolia_Neolithic, Anatolia_Chalcolithic, Levant_Neolithic, Levant_Chalcolithic, Zagros_Neolithic, Zagros_Chalcolithic, Mesopotamia_Neolithic, Caucasus_Neolithic and Caucasus_Chalcolithic) as well as two Neolithic Europeans: Spain_EN and Greece_Neolithic, referred to as the full qpAdm model. This analysis was conducted over 474,731–578,969 SNPs. Third Intermediate Period ancestry modelling We estimated the ancestry proportions of the two Third Intermediate Period Egyptians20 using North African and West Asian populations who lived between the Old Kingdom and Third Intermediate periods as potential sources (NUE001, Anatolia_BA, Levant_BA, Iran_BA and Caucasus_BA), as well as a Bronze Age Greek population (Greece_Minoan). We also added Morocco_MN and Mesopotamia_N to test whether the Third Intermediate Period Egyptians share a closer ancestry with NUE001 or a source more related to one of these two ancestries present in NUE001. This analysis was conducted on 290,262 SNPs. Present-day Egyptian genome ancestry modelling We directly estimated the proportion of NUE001 ancestry in present-day Egyptians164 as a whole or each individual separately, as well as ancestries from North African (Morocco_MN), West Asian (Caucasus_BA, Iran_BA and Levant_BA) and European (Greece_Minoan) populations, as well as East and West Africa (Ethiopia_4500BP.SG and Congo_Kindoki_Protohistoric). For each region, we selected the representatives closest to NUE001’s lifetime. Anatolia_BA was removed from the list of West Asian groups because its inclusion led to a substantial drop-off of genomes having at least one model passing P = 0.05. This analysis was conducted on 767,305 SNPs. Ancient East African ancestry modelling We estimated the ancestry proportion in ancient East African43,54,55,144 using both NUE001 and Levantine Chalcolithic genomes as competing sources for the Eurasian-like component. We used as potential left sources NUE001, Levant_Chalcolithic, Ethiopia_4500BP.SG, Dinka.DG, Congo_Kindoki_Protohistoric and South_Africa_2200BP.SG. For this model, the following fixed right groups were used: Chimp.REF, Latvia_HG_UDG, USA_Ancient_Beringian.SG, Vanuatu_400BP_UDG, Japan_HG_Jomon_UDG and China_NEastAsia_Coastal_EN_UDG. This analysis was conducted on 141,323–350,110 SNPs. f 4 statistics f 4 statistics of the form f 4 (NUE001, Morocco_MN; X, Ju_hoan_North.DG) and f 4 (NUE001, Morocco_MN; Mesopotamia_N, X), X being the non-North African groups used in the full qpAdm model, Palaeolithic Levant or Palaeolithic Anatolia, were estimated to confirm the probable source of admixture in the Nuwayrat genome when compared with the Middle Neolithic Moroccan group. f 4 statistics was computed using ADMIXTOOLS 2 (ref. 42) with the option maxmiss = 0.1. We restricted the analysis to genomes with both transitions and transversions (half/plus UDG-treated libraries or single-stranded libraries called using SequenceTools115 --singleStrandMode). This analysis was conducted on the 1240k dataset on 280,544 SNPs. Imputation The genotypes of the Nuwayrat genome were imputed together with 200 ancient genomes from North Africa and West Asia associated with the Palaeolithic, Neolithic and Bronze Age culture (Supplementary Table 3). We restricted the imputation to whole-genome sequencing data greater than 0.5× coverage or 1240k SNP capture data greater than 2× coverage, following recommendations from Sousa da Mota et al.172. First, we called genotypes using bcftools v.1.19 (ref. 173) with the commands bcftools mpileup with parameters -I -E -a ‘FORMAT/DP’ --ignore-RG and bcftools call -Aim -C alleles. We then imputed the missing genotypes using Glimpse v.1.1.0 (ref. 174). First, we used GLIMPSE_chunk to split chromosomes into chunks of 2 Mb and with a 200-kb buffer region. Second, imputation was performed with GLIMPSE_phase on the chunks with default parameters --burn 10, --main 10 and --pbwt-depth 2, with 1000 Genomes175 as the reference panel. We then ligated the imputed chunks with GLIMPSE_ligate. To remove transitions caused by post-mortem damage before imputation, for the genomes generated with UDG treatment, we first hard-trimmed the first and last three base pairs of each read and removed CpG sites, and for the genome generated without UDG treatment, we removed all transition sites after SNP calling. We finally restricted the imputed genotypes to those with genotype probability ≥0.99 and minor allele frequency ≥0.01 using the command bcftools filter -i ‘MAX(FORMAT/GP)>=0.99 && INFO/RAF>=0.01&&INFO/RAF<=0.99’ --set-GTs ‘.’. The imputed dataset was used for phenotype prediction (see below) and admixture dating using DATES (Supplementary Information section 4 and Supplementary Table 11). Phenotype prediction The genotypes responsible for skin, hair and eye colour prediction were investigated using the HIrisPlex-S system176,177,178 using the imputed genotypes. Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.