Ancient Maltese genomes and the genetic geography of Neolithic Europe

Summary Archaeological consideration of maritime connectivity has ranged from a biogeographical perspective that considers the sea as a barrier to a view of seaways as ancient highways that facilitate exchange. Our results illustrate the former. We report three Late Neolithic human genomes from the Mediterranean island of Malta that are markedly enriched for runs of homozygosity, indicating inbreeding in their ancestry and an effective population size of only hundreds, a striking illustration of maritime isolation in this agricultural society. In the Late Neolithic, communities across mainland Europe experienced a resurgence of hunter-gatherer ancestry, pointing toward the persistence of different ancestral strands that subsequently admixed. This is absent in the Maltese genomes, giving a further indication of their genomic insularity. Imputation of genome-wide genotypes in our new and 258 published ancient individuals allowed shared identity-by-descent segment analysis, giving a fine-grained genetic geography of Neolithic Europe. This highlights the differentiating effects of seafaring Mediterranean expansion and also island colonization, including that of Ireland, Britain, and Orkney. These maritime effects contrast profoundly with a lack of migratory barriers in the establishment of Central European farming populations from Anatolia and the Balkans.


In brief
Using high-coverage genomes from ancient Maltese, Ariano et al. show that they had signals of both recent and more distant inbreeding, indicating small population size. Using diploid genomewide data, the authors explore the genetic geography of Neolithic Europe which, echoing modern patterns, is shaped by physical topography, especially its seascapes.

INTRODUCTION
The importance of sea travel in prehistory is clear from the rapid westward spread of agriculture from its origins in the Near East along the Mediterranean littoral, including its very early appearance in Cyprus circa (c.) 10,600 years ago. 1 However, the consideration of seascapes in prehistory has varied, with a biogeographical view emphasizing the sea as a barrier and, alternately, a view that posits seaways as efficient corridors of connectivity. 2,3 Ancient genomics has confirmed the demic, or migratory, nature of Neolithic expansion but has also given some illustrations of retardation of seaborne genetic exchange. For example, the Sardinian Bronze Age population was unaffected by an influx of Steppe ancestry that changed the genomes of contemporaneous mainland Europeans, 4,5 and Irish Mesolithic genomes show the signatures of small population size, which were absent in corresponding continental hunter gatherers (HGs). 6 The first settlements in the Maltese islands were Neolithic, dated from the sixth millennium BC. These developed through a series of cultural phases, with some material indications of external connectivity, 7 but faded from 3600 BC when pottery and architecture started to show distinctive features. 8 One example was the development of multi-chambered rock-cut tombs, such as that at Xag-ra circle, Gozo ( Figure 1). 7 This monumentalized underground tomb yielded the remains of hundreds of individuals 9 and underwent remodeling and enlargement until around 2500 BC when it was abandoned, possibly as part of a wider population decline or replacement.
To examine the demography of Late Neolithic Malta, we sequenced genomes from Xag-ra circle. The elucidation of fine structure among closely related groups such as European Neolithic populations is challenging and requires the resolution afforded by genealogical methods. 10 Therefore, to examine these in a wider context, we additionally imputed genomewide diploid genotypes from published ancient genomes and assessed haplotype sharing within and between genomes to estimate genetic geography and demographies across Neolithic Europe.

RESULTS AND DISCUSSION
Genomes from a south Mediterranean island The retrieval of ancient genomes from warm climates is highly challenging, and the island of Gozo in the Maltese archipelago is one of the southernmost contexts in Europe ( Figure 1). However, from nine human petrous bone and tooth samples from the Late Neolithic Xag-ra Circle excavation, three yielded excellent endogenous DNA content (13%, 17%, and 39%; Table 1). This likely reflects enhanced preservation within this underground limestone cave burial complex (hypogeum), and these samples, Xaghra5, 6, and 9, were shotgun sequenced to an average genome-wide coverage of 1.243, 0.983, and 7.523, respectively.
Malta was one of the final regions of Europe to be inhabited, with little evidence of human presence prior to the arrival of Neolithic communities, which were established on the archipelago by 5500 cal. BC. 11 These were associated with a developed style of impressed pottery (G-ar Dalam ware) that represented a regional variant of Sicilian and southwestern Italian ceramics. Article Accordingly, we find the genomes from Xag-ra Circle share highest levels of drift with the Early Neolithic populations of Italy and Greece, followed by Middle Neolithic and Chalcolithic populations from Italy and Sicily, as estimated using outgroup f 3 -statistics (Data S2F; Figure S2). Levels of Western HG (WHG) admixture have been shown to vary across European Neolithic samples, [12][13][14][15][16][17] particularly through time. To examine levels of WHG ancestry within our Neolithic sample, we applied the qpAdm method to each site, binning genomes into 500-year intervals. We observe WHG ancestral components to increase significantly with time ( Figure S3; r 2 = À0.52, p value = 2.8eÀ4). Interestingly, the Xag-ra Circle site shows the lowest amount of HG ancestry (6.8% ± 2.5%) among other groups from the Later Neolithic ( Figure S3). This may reflect a shielding by its island context from the dissemination of admixtures with persisting WHG populations that widely influenced mainland populations and which have been estimated to occur as late as 3800 BC. 17 This resonates with observations from Sardinian populations, which show a constant degree of WHG ancestry stretching through the Neolithic to Bronze Age periods. 4,5 Using D-statistics, we also tested for gene flow related to North African, Caucasus HG, Neolithic Iranian farmer, and Yamnaya-steppe groups into the Maltese populations, to the exclusion of the Greek and Italian Early Neolithic. We obtained no statistical evidence for admixture (Data S2). The imputation of diploid genotypes from low coverage shotgun sequence data has been successfully utilized for the characterization of fine-scale structure and patterns of inbreeding in ancient populations. 6,18,19 We applied an imputation pipeline (STAR Methods) to the Maltese samples using Impute2, 20,21 as well as to 117 individuals for which sufficient shotgun sequence (>0.43) was available (Data S1A). The resulting diploid genotypes were merged with relevant ancient Italian samples from a published imputed dataset. 22 We also extended an imputation pipeline to individuals that had been sampled with a SNP capture protocol 13,14,23 using Beagle v.4.1, 24,25 achieving an accuracy in the predictions of heterozygous genotypes of 95% ( Figure S6). After excluding 4 samples with high numbers of missing genotypes, this gave a final comparative dataset of 271 Neolithic and 86 HG ancient individuals from western Eurasia 5,6,12,13,15,16,18,19,22, (Data S1A). Comparisons of runs of homozygosity (ROHs) estimates, using diploid highcoverage data and the alternate imputation pipelines for individuals with both shotgun sequence and SNP capture data, show very high concordance (r 2 = 0.99; Figure S4) and validate the combined analysis of this dataset. Comparisons of identity-bydescent (IBD) scores between the data types also show no evidence of bias ( Figure S4).
Xag-ra circle genomes show outlying homozygosity levels and a historically restricted population size Genome-wide diploid data allow haplotype-based assessments of population diversity-specifically, the distribution of shared ancestry within genomes, using ROH, and the distribution between individuals by identifying shared tracts that are identical by descent. ROH analysis shows outlying behavior by the Maltese genomes. Xaghra9 has the second most extreme levels of long ROH (>5 cM) yet reported in prehistory: an assertion secured by its high genome coverage ( Figure S5) and a confirmatory analysis using a second-analysis method (using RO-Han 52 ), which estimated 19.12% of the genome under ROH. This is only exceeded within an individual deposited in an Irish passage tomb (NG10, Newgrange10), who was the offspring of a first-order consanguineous union. 6 However, Xaghra9 has a ROH size spectrum that has less skew toward very long tracts of identity (>15 cM; Figure 2A).

Article
To explore this signal, a range of consanguineous parentages were simulated, and the number of ROH segments with the total fraction of the genome in these ROH (F ROH ) were plotted and compared with ancient individuals ( Figure 2B). Unlike NG10, Xaghra9 falls at the edge of the distribution seen for matings between first-degree relatives and may result from a more complex combination of multiple inbreeding loops within his genealogy. However, this is similar to Israeli Chalcolithic sample I1178 43 (F ROH = 0.16; Data S1D), who was previously identified in a different analysis as a possible product of brother-sister or parent-offspring consanguinity. 53 Consequently, we do not assert a precise scenario for the parentage of Xaghra9. To focus on very recent inbreeding, we repeated this analysis twice, considering only ROH segments longer than 10 cM and then longer than 15 cM; for each, Xaghra9 remains at the edge or outside the sibling mating cluster. Given the small size and relative isolation of Gozo island, it is possible that the inbreeding loops that gave rise to the Xaghra9 genome are the result of both recent genealogical inbreeding and a historically small ancestral population size. This interpretation is supported by the observation of less pronounced but relatively inflated levels of the fraction of the genome in ROH in the other two Maltese genomes (Xaghra5, Xaghra6; Figures 2A, 2B, and S5A), one of which predates Xaghra9 by $400 years. The values for these two samples are more typical of those found in European HGs, who maintained smaller population sizes than later farming populations ( Figure 2B). To investigate further, we used levels of ROH within a range of 4-20 cM and a maximum likelihood framework 53,54 to estimate effective population size, giving a total of 515 (95% confidence interval [CI] 397-633) individuals.
We also calculated effective population size for the Xag-ra population using the software IBDNe, 55 which leverages patterns of IBD sharing between individuals. For comparison, we included other European Neolithic sites with more than 90 IBD segments shared between individuals in total. Xagra, and to a lesser extent the remains from the Tomb of the Eagles on Isbister in the Orkney islands, show recent dips in population size, with the Late Neolithic Maltese sample giving a 30-generation average of only 382 individuals ( Figure 3A).
Thus, these preserved Maltese samples show a genomic signature of an unusually small and restricted population, a signal which is distributed over a period of at least 400 years. Interestingly, the two later individuals (Xaghra5 and Xaghra9) derive from a turning point in Maltese prehistory c. 2450 BC, with a reducing density of radiocarbon dates 56 and marked worsening in diet and nutritional status. 57 A long-term trend toward increasing aridity and thinning soils that began as early as 5500 BC 58 seems to be driving these changes, implying the Late Neolithic population was less than the Early Neolithic carrying capacity estimate of two or three thousand individuals for Gozo island (67 sq km). 58 This is only a small multiple of our calculated effective population size values, which are therefore See also Figure S5 and Data S1D. Article not surprising. However, these estimates suggest isolation, with mating networks largely confined within the island's shores. Several strands of evidence suggest our sample is representative of the wider Neolithic community on Gozo. First, the age profile of Xag-ra burials coincides closely with expectations of the mortality rates of a full early farming community, namely high infant and adolescent mortality and a relatively equal balance of adult males and females. 59 Second, the spatial analysis of the mortuary remains suggests a rich and elaborate treatment of the burials as one community, 9,60,61 -12166603600450-2819404725035 and finally, the chosen samples are drawn from different parts of the site and span the entirety of its use. Archaeological evidence for overseas communication with Malta in this period is mixed. Some products, such as obsidian, types of chert, and polished stone were definitely imported. 60,62 However, these tend to be small, of high prestige value, and have a finished state when they appear, suggesting they may not have been accompanied by a substantial volume of human traffic. Moreover, the means of cultivation of crops, raising of animals, and construction were local in nature, consistent with a degree of insularity. See also Data S1B.

Article Body mass index analysis of Neolithic populations
The carving and circulation of apparently obese human figurines was a marked feature of the late Maltese Neolithic, 8 perhaps mirroring an unusual genetic predisposition within a restricted gene pool. Accordingly, we performed a polygenic risk score analysis on body mass index using the summary statistics from the UK Biobank dataset, but found that the three Maltese Neolithic individuals sampled do not give atypical risk values compared with other Neolithic individuals ( Figure S8).
Haplotype sharing within and between Neolithic sites suggests restricted island population sizes and seaborne founder effects Shared IBD is sensitive to recent common ancestry and, because it is a genealogical rather than a frequency-based method, 63,64 it may be less skewed by factors such as the differences in levels of HG ancestry that are known among European Neolithic populations. [12][13][14][15][16] Figure 3C (Data S1B) shows a heatmap of the average IBD length (R2 cM) observed between and within European Neolithic archaeological sites with more than one imputed genome, after filtering for related individuals.
We observe the highest within-site values for samples from small islands, with Xag-ra (Malta) producing the most extreme result, followed by Ansarve (Gotland), Holm of Papa (Orkney), and Isbister (Orkney), supporting restricted population histories for insular Neolithic societies. Figure 3B plots the averaged values for different geographical regions and reveals an additional trend of higher within-group IBD sharing in the north and west of the continent relative to the south and east. This geographical difference also manifests in patterns of between-site sharing ( Figure 3C), with three distinct regional clusters apparent. The Basque region, situated between the Atlantic Ocean and the western Pyrenees mountains, shows extremely inflated values between Later Neolithic sites, implying a degree of geographic isolation. Close genealogical ties are also seen across Britain and Ireland, consistent with a seaborne colonization of the islands derived from a single or closely related founder populations. Finally, we observe French sites clustering together, within which extreme sharing is observed between two Early Neolithic sites from Southern France, potentially reflective of the enclave colonization process that characterized Neolithic expansion across the Mediterranean. To explore this signal further, we considered three sites from the earliest horizon of the Spanish Neolithic (c. 5500-5000 BC), previously excluded given only a single sample was available from each. Surprisingly, despite the large geographic distances between them (Data S1A), these three individuals show very high levels of sharing with one another and with the Mediterranean French sites, despite large differences in their HG ancestral contribution (Figure S3). This implies a population size restriction accompanied Neolithic migration into the Western Mediterranean.

Neolithic genes mirror geography
To explore the potential impact of maritime colonization and continental topography on Neolithic genetic structure, we carried out principal component analysis (PCA) on a matrix of pairwise IBD sharing between individual imputed ancient individuals (Figure 4), as well as ChromoPainter and clustering using fineS-TRUCTURE analysis 65 ( Figure 5). In addition to these haplotype-based methods, we also applied an allele frequency-based approach (EEMS, estimated effective migration surface 66 ).
Results from each show a convergence on the existence of three clusters: first, the Western Mediterranean, including Iberian, French, and Sardinian individuals; second, the Eastern Mediterranean, featuring Greek, Balkan, and Anatolian individuals as well as Central Europeans; and third, the British and Irish archipelago. These are visible as blocks in the IBD heatmap of Figure 3C and form three apices of variation in the PCA (Figure 4). They also form separate primary branches in a fineSTRUCTURE tree ( Figure 5). 65 Intermediate samples are also intermediate in geography. For example, in the PCA plot, which visibly mirrors geography (Figure 4), Northern French samples are placed close to Iberians but also stretch toward the British and Irish cluster. Also, the mid-Mediterranean samples from Sardinia, Malta, Sicily and Italy fall between the western and eastern poles.
Neolithic populations migrated through Europe via two major routes, an overland transfer into Central Europe and a maritime dissemination along the Mediterranean coast. 67 The most striking feature in our analyses are the contrasting outcomes of these two processes. Particularly, there is minimal distinction between central European individuals and their source populations in the Balkans and Anatolia, whereas the separation of western European individuals from those in the southeast forms the primary divide in the data. This supports a model of agricultural expansion into Central Europe from the Balkans that involved substantial numbers of migrants and strong backward communication during the dissemination of the Linearbandkeramik (LBK) complex, with populations remaining relatively well connected throughout the Neolithic period.
To explore further, we also EEMSs using a stepping-stone model and a distance matrix computed from allele frequencies 66 ; Figure 6 shows cold-and hotspots of estimated migration rates within Neolithic Europe. The communication corridor between Anatolia, the Balkans, and Central Europe is the most striking feature of this analysis and contrasts strongly with east-west barriers in the Mediterranean sea, the Alpine region, and further north where the two Neolithic migratory streams are purported to meet. 51 In common with the other approaches, EEMS does not take account of temporal differences among samples, which would be expected to be a differentiating factor. For example, the barrier between English samples and the continent might be less pronounced with the addition of more contemporaneous French genomes. However, we assert that the major divisions are explained at least partially by geography. These correspond with those that emerge in the haplotype-informed fineSTRUCTURE analysis, where sample dates are also plotted ( Figure 5). From this it is clear that genomes separate into different groups despite overlapping contemporaneity across the basal branches. Also, there are considerable temporal differences within clusters, particularly among the samples in the Anatolian-Central European high communication corridor.
The rapid Neolithic colonization of the Western Mediterranean from the east was associated with the impressed cardial complex and likely took place through iterative coastal placements along the northern maritime littoral. 68 Article is required to explain the speed of agricultural spread, which was significantly faster than that seen in Central Europe. 67,70 Our results accord with a limited capacity of sea craft used in this cabotage, which likely restricted pioneer numbers and subsequent backward exchange. We infer that the observed east-west genomic distinction derives at least partially from this foundational process, as earlier individuals plot toward the extremes in the PCA graph with mid and Late Neolithic individuals showing a more central tendency ( Figure 4C). The sharp divide between eastern and western Europe echoes the analysis of French and neighboring Neolithic genomes by Rivollat et al., 51 who also identify that the two Neolithic streams differed in their degree of ancestral admixture with European HGs. However, this difference in ancestry is less marked in comparison with earlier western genomes, for example, those of the Iberian Early Neolithic ( Figure S3).
British and Irish populations form a sister grouping to the Mediterranean Neolithic in the second fineSTRUCTURE branching ( Figure 5) and visibly show IBD affinity ( Figure 3C), according with prior assertions that they primarily owe their origins to this southern migratory stream. 6,15,41,49 However, their maritime separation is mirrored by a degree of cluster distinction (Figures 3  and 4) and an estimated migration barrier ( Figure 6). Interestingly, Irish and mainland British individuals do not separate from each other as clusters in any of our analyses, supporting shared elements of a rapid foundation process c. 3800 BC. 71 This is an additional indication of the absence of significant batch effects, as the British were imputed from SNP capture data and the Irish from shotgun-sequenced libraries. However, fineS-TRUCTURE confirms the emerging distinctiveness of (SNPcaptured) Orcadian individuals, as well as that of Basque Late Neolithic sites ( Figure 5), also captured in patterns of IBD sharing ( Figures 3C and 4). An additional marker of separation is that Orkney islander ancient genomes have also recently been found to show unusual majority retention of male lineages across the Neolithic-Bronze Age transition, 72 a feature unique within Northern and Central Europe.

Conclusions
Basque, Orcadian, and Irish distinctiveness emerged in pioneering studies of modern human genetic variation, [73][74][75][76] and genome-scale investigation has compellingly recapitulated the geography of Europe in PCA, particularly its maritime features. 77 It is striking that these same features emerge independently within data from an earlier genomic era in the same continent, speaking to the repeated shaping of genetic variation by the same physical topography, particularly its seascapes. One of the great debates of prehistory is the level of maritime

OPEN ACCESS
Article connectivity during the course of millennia and how that connectivity interacted with marine technology and cultural response. 2 We suggest that relationships among ancient European populations indicate that sea travel was one driver of genomic differentiation during the establishment of the Neolithic. On a wide scale, multiple analyses highlight the genetic separation between Western Mediterranean sites and their source Eastern

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:  Computed using a stepping-stone model and imputed allele frequency data, migration rates are plotted as log 10

DECLARATION OF INTERESTS
The authors declare no competing interests.

INCLUSION AND DIVERSITY
We worked to ensure diversity in experimental samples through the selection of the genomic datasets. One or more of the authors of this paper self-identifies as a member of the LGBTQ+ community. One or more of the authors of this paper self-identifies as living with a disability. One or more of the authors of this paper received support from a program designed to increase minority representation in science. While citing references scientifically relevant for this work, we also actively worked to promote gender balance in our reference list.

ll
OPEN ACCESS cleaned using a drill bit. A triangular wedge section of the otic capsule region of each petrous bone and the root of each tooth were cut using a Dremel diamond wheel saw. Each sampled bone part was pulverised in a Mixer Mill MM 400 (Retsch). An aliquot of $0.1g of this bone powder was used for DNA extraction, and the rest of the powder was stored in a separate tube. The DNA extraction procedure followed the same protocol described in Yang et al. 95 with modifications presented elsewhere. 96 One sample subsequently sequenced at high coverage was re-extracted using an initial washing step by 0.5% bleach solution as described in Boessenkool et al.. 97 Radiocarbon analysis Date ranges have been estimated using the 95% confidence interval of Bayesian chronological models of 117 radiocarbon dates from the site and their stratigraphic relationships. 7 The three sequenced samples all derive from the megalithic burial Circle on the Xag-ra plateau between the temples of _ Ggantija and Santa Verna, excavated between 1993 and 1994. The oldest sample (Xaghra6) derives from a deeper stratified area of stacked burials that also contained rich ceremonial objects higher in the stratigraphy. A burial in the same layer as Xaghra6 was radiocarbon dated to 2900-2650 BC (OxA-27837, 4198±26 BP). 7 The two later samples were found in shallower deposits to the west. Xaghra5 was part of a display area of initially articulated human remains placed with portable figurines that was intentionally dismembered, most probably, at least in part to the slightly deeper location of Xaghra9 slightly to the north. Xaghra6 was placed as the use of the site began to intensify whereas the other two samples date to the period of most intensive activity some four to five hundred years later, with 23 radiocarbon measurements from material associated with these samples spanning approximately 2550 to 2350 BC.

Library preparation
The initial screening of each sample and blank controls was performed by constructing a double-stranded DNA NGS library, priorly treated with Uracil-DNA-glycosylase (UDG), using the method outlined in Meyer and Kircher 98 and modified as described in Gamba et al. 12 Libraries were amplified with AccuPrime Pfx Supermix (Life Technology) using 12-14 cycles of PCR, assigned with unique indexes and quantified with a TapeStation 2200 (Agilent Technologies). The same libraries were also used for further amplifications required for high coverage sequencing.

DNA sequencing
The initial screening to assess the endogenous DNA was performed by sequencing all the libraries with the Illumina HiSeq 2500 platform (100bp SE) at Macrogen (Republic of Korea). Subsequently, 3 samples with high endogenous DNA were further sequenced to high coverage using the HiSeq 2500 Illumina platform (100bp SE) at Macrogen (Republic of Korea). One sample was further sequenced using NovaSeq (50bp PE) Illumina platforms at TrinSeq (Ireland).

Reads processing
For samples sequenced in single-end mode, reads were trimmed of their adapters and filtered based on their length using the software cutadapt v.1.9.1 79 (cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -O 1 -m 34). For paired-end libraries, adapters were trimmed and reads were filtered using AdapterRemoval v2.1.1 80 (-trimns -trimqualities -minquality 25 -collapse). Reads that passed these qualities and length filters were aligned to the human reference genome (hg19/GRCh37) with the mitochondrial sequence replaced by the Revised Cambridge Reference Sequence (rCRS, NC_012920.1) using the software BWA v.0.7.5a 81 with relaxed parameters (-l 1024 -n 0.01 -o 2). Aligned reads that came from PCR duplication or with a mapping quality below 20 were removed using the software SAMtools v.1.7 82 and Picard Tools v.1.101 (http://broadinstitute.github.io/picard/). The coverage of each completed aligned file was calculated using the tool Qualimap v.2.1.1. 83 Indels were locally realigned using The RealignerTargetCreator and IndelRealigner tools from GATK v.2.4. 85 Additionally 2bp were soft clipped at the start and end of each read.

Contamination estimation and sex determination
To determine the sex of each sample we applied two methods, one outlined in 84 and the other described in Cassidy et al. 6 In both methods, the amount of reads aligned on the X chromosome versus the autosomal genome was used to estimate the sex of an individual together with a confidence interval. We only considered sex assignments where both methods agreed. For three Maltese samples analysed in this study we estimated contamination using the haploid information contained in the mitochondrial genome and in the X chromosome for two males, applying the same method outlined in. 15 Population structure analysis Pseudohaploid genotypes were called at approximately 600,000 autosomal sites from the Human Origins panel 27 for the same set of ancient samples used in the IBD analyses plus 19 other ancient samples representative of hunter-gatherers, Bronze Age and Neolithic farmers populations. 13,14,18,99 Read bases were determined at each site using the Pileup tool from GATK v2.4, 85 filtered for a quality of 30, with bases not matching either the reference or alternate allele removed. A single base was then randomly selected to generate the pseudohaploid genotype. This ancient dataset was then merged with a subset of the Human Origins panel from Article Western Eurasia using the software PLINK v1.9. 86 A Principal Component Analysis (PCA) was then carried out on the 604 modern individuals from Human Origins, with the genetic variation of the ancient samples projected onto this using the SmartPCA v.16000 algorithm implemented in EIGENSOFT 87,100 with parameters (killr2: YES, r2thresh: 0.2, numoutlieriter: 0, lsqproject: YES, autoshrink: YES) ( Figure S1).

F-statistics
Using the same set of ancient samples described in the previous paragraph and transversion sites only from the ''1240K'' panel, 13,14,23 we estimated the amount of drift that the Maltese shared with each other population using the outgroup-f3 statistics 88,101 method implemented in the ADMIXTOOLS package v.7.0.2. 88 This analysis was carried out in the form of f3(Mbuti; Ancient Maltese, X) where X represents different populations tested (Data S2F; Figure S2). The outgroup population, Mbuti, is represented by four individuals collected from the SGDP dataset. 102 To test for Admixture with African populations we used D statistics. 103 Four ancient North African representatives were selected from Fregel et al. 104 (Data S2A). Tests were constructed in the form of: D(Chimp, Ancient North-Africa, Malta Late Neolithic, X) where X represents Neolithic populations that fall close in the PCA to the ancient Maltese ( Figure S1; Data S2D).
Similarly to test for admixture between the Maltese and Caucasus hunter-gatherer (CHG) or Steppe populations we built our D statistics test in the form of D(Mbuti, CHG/Yamnaya, Malta Late Neolithic, X). In this analysis the CHG population is represented by two individuals published in Jones et al. 18 (Data S2B and S2C).
To estimate the amount of WHG ancestry we used the method qpAdm Haak et al. 14 We first divided the individuals into groups according to their archaeological site of origin. Each group was furthermore subdivided in bins of 1000 years and only sub-groups with at least 2 individuals were considered for this analysis. The reference group was comprised of the following genomes: (Mbuti.DG, Ust-Ishim, MA1, Villabruna, GoyetQ116-1, Han.DG, Papuan.DG, Mixe.DG, Karitiana.DG, AHG, Iran_Neolithic, CHG, EHG). The source population are Anatolian_Neolithic represented by individuals from Barcin and WHG individuals represented by Loschbour and KO1 (Data S2E). Only groups with a p-value higher than 0.05 were included in Figure S3.

Genotype imputation
From samples that had been screened using an in-solution target capture method we selected 231 published genomes for imputation with a reported coverage on target regions of at least 0.6X and 650K SNPs called from the 1240K panel. To increase the number of samples from Neolithic Sardinia we also included 4 samples with a coverage higher than 0.6X and at least 460K SNPs safely called from the 1240K panel. Before imputation we selected a set of approximately 6.2 million SNPs to be called on our target dataset using the 1000 Genomes Project (1000G) resource as reference, filtered for individuals of African origin (defined with the AFR label) and with a minor allele frequency of 5%. Variants were called using the tool UnifiedGenotyper in GATK v2.4 85 program with parameters (-output_mode EMIT_ALL_SITES, -genotyping_mode GENOTYPE_GIVEN_ALLELES). The VCF files created were then split first by chromosome and then by windows of 1 Mb. Genotype imputation was performed on approximately 28 million variants using the tool Beagle v.4.1 24,25 with a reference dataset of 1843 modern individuals of non-African origin from the 1000 Genomes project. The program was run in multi-thread mode taking advantage of the Irish Centre for High-End Computing (ICHEC) cluster. The genetic map used was taken from the Beagle website (http://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/). The imputed VCF files were filtered for SNPs only and genotype probability of 0.95 using bcftools v.1.3 82 and PLINK v1.9 (-vcf-min-gp 0.95) 86 obtaining 25.8 million variants.
After completion of imputation, four samples with high genotype missingness (>=0.1) were removed from subsequent analyses. Separately we also selected 120 WGS samples with a coverage of at least 0.4X to impute using the Software Impute2. 20,21 For these samples, and similarly to the SNP capture imputation, variants were called from the 1000 Genomes project 78 using the tool UnifiedGenotyper in GATK v.2.4 85 using the same parameters. The whole 1000 Genomes project dataset was used as reference for the genotype imputation. Prior to imputation transition SNPs were excluded from this dataset resulting in calling of approximately 28 million. The VCF file was then split first by chromosome and then in windows containing 15000 markers. For each input file, the program Impute2 was called using the parameters ( -Ne 20000 -buffer 500 -allow_large_regions -k 400 -k_hap 2000). After imputation we filtered for genotype probability higher than 0.99 (GP > 0.99) resulting in 77.8 million SNPs.
Finally this combined resource was merged with 21 low coverage shotgun sequenced genomes which had been previously imputed. 22 Between WGS and SNP capture samples we obtained a final resource of 357 unique imputed diploid genomes (Data S1A; Figure S6).
After merging these resources we then tested for differences in genotype missingness between datasets. To do so we first considered a set of 12 million SNPs common across all three datasets. We then calculated the missingness for each dataset and averaged across samples. We observed a genotype missing for the SNP capture and WGS imputed respectively of 12% and 13.5%.

ROH and inbreeding analysis
To estimate the inbreeding coefficients of our imputed samples, we used a measure based on the proportion of the genome that is homozygous-by-descent (runs of homozygosity that are identical by descent), as employed in Cassidy et al., 6 and labelled here as F ROH . Separately, the hunter-gatherer and Neolithic farmer datasets were filtered for genotypes missingness and minor allele frequency using PLINK v1.9 (-geno 0.02, -maf 0.05, -indep 50 2 2) obtaining respectively 51,289 SNPs and 41,426 SNPs. Using this set of SNPs we then identified ROH segments using PLINK v1.9. with the same parameters used in Gazal  Article (-homozyg-window-het 0 -homozyg-snp 50 -homozyg-kb 1 -homozyg-density 5000 -homozyg-gap 5000). Physical measures were converted to centiMorgans (cM). The total length of the genome in ROH above this threshold divided by the length of the autosomal part was used to estimate the F ROH coefficients. 106 To assess the concordance between samples imputed from different sources we compared F ROH estimates obtained from imputed SNP capture with those calculated using imputed data from WGS data available for the same samples. We considered the same set of SNPs in both data types. For two samples that were whole genome screened and where the coverage was sufficiently high we also estimated the F ROH coefficients using diploid genotype calls. For these two samples we applied the same protocol described in the imputation accuracy paragraph. In brief diploid genotypes with a depth below 10 or higher than 30 and a quality below 30 were excluded. As shown in Figure S4 there is no visible deviation of substance between the measures.
For the sample, Xaghra9, which has sufficient coverage, we ran the software ROHan 52 to validate our inbreeding results. As suggested by the software we first used the program bam2prof with different threshold values ( -length 5, 10, 15, 20) to account for postmortem deamination damage. We then run the program rohan using transversion only (-tvsonly).

Pedigree simulation
To better understand the degree of relatedness between the parents of inbred samples we simulated different pedigree scenarios using dummy genotypes. We started from the same dataset described in the previous paragraph and we filtered for genotype missingness and minor allele frequency. This filtered resource was split by chromosome and then re-phased using SHAPEIT v.2.r837. 90 After phasing we filtered for linkage disequilibrium with plink using (-indep 50 2 2) and selected a common set of SNPs. 21 Irish imputed individuals published in Cassidy et al. 6 were selected from this dataset as founders to build the simulated pedigrees. This set of founders were not influenced by inbreeding, relatedness, population structure, or recent change in population size. This dataset was then used as input for the software ped-sim 91,92 with a refined genetic map taken from Bh erer et al. 107 Three different inbreeding scenarios were tested: -First degree: siblings and parent-offspring -Second degree: uncle-niece/aunt-nephew and grandparent-grandchild -Third degree: first cousins and great aunt-great nephew/great uncle-great niece Each of these scenarios was simulated 400 times using random sampled founders. ROH segments were found using PLINK with the same parameters described in the previous section (-homozyg-window-het 0 -homozyg-snp 50 -homozyg-kb 1 -homozyg-density 5000 -homozyg-gap 5000) and inbreeding coefficients estimates were also obtained using the same pipeline for both simulated and real genomes.

IBD analysis
In this work, we applied the software IBDSeq vr1206 94 to the unphased dataset to identify segments of the genome inherited by recent common ancestors (identical by descent) in European Neolithic samples. Genotype missingness and minor allele frequency filters were applied to the imputed dataset using the software PLINK v.1.9 (-geno 0.02, -maf 0.05). Related individuals with a relatedness estimated by the software KING v.2.2.6 105 higher than 4th degree relatives were also removed from analyses obtaining 258 unrelated samples. Filtered files in PLINK format were converted to VCF using the option (-vcf) in PLINK v1.9. and used as input to the program IBDSeq with parameters (errormax=0.005 and LOD >= 3; 108 ). IBD segments shorter than 2 cM were excluded following the advice of Browning and Browning. 94 To test that no systematic bias was present between types of data, we compared the results obtained from those samples where it was possible to impute genome wide calls using both WGS and SNP capture data. We used a common set of SNPs for both data types that were pruned for genotype missingness and minor allele frequency, obtaining approximately 900 thousands markers per comparison. This set of common SNPs was then used to calculate the total amount of IBD that each sample type, WGS or SNP capture, shared with the rest of the dataset. As shown in Figure S4 correlation and variation around the 1:1 plot line indicate no systematic bias between captured and WGS imputed data Population size estimates To estimate the effective population size we used the IBD information obtained from IBDSeq as an input for the software IBDNe v.23Apr20.ae9. 55 This software was run for 50 generations with default settings and only for groups that shared at least 90 IBD segments longer than 2cM. An estimate of population size for each group was calculated by taking the harmonic mean over 25 generations (from 5 to 30).
Separately we also estimated the effective population size of our Maltese group using the software hapROH v0.3a4. 53 First we excluded the highly inbred sample Xaghra9 from this group. For the remaining two imputed samples(Xaghra5 and Xaghra6), diploid genotypes were downsampled to ''1240K'' SNPs panel and ROH were called with plink similar to what is described above (-homozyg-window-het 0 -homozyg-snp 50 -homozyg-kb 1 -homozyg-density 5000 -homozyg-gap 5000). For each of the two Maltese samples the ROH results were then used to estimate the effective population size using the function ''MLE_ROH_Ne'' from the hapROH package using the parameters (min_len=4, max_len=20, ne=10000, bin_range=[0.04, 0.5], nbins=1000, error_model=False). Article Chromopainter/fineSTRUCTURE To investigate fine-scale population structure in our imputed dataset we used the software fineSTRUCTURE v2. 65 The same set of unrelated samples used in the IBDseq analysis were used for this analysis. These ancient imputed samples were filtered for genotype missingness and minor allele frequency using the software PLINK v.1.9. with parameters (-geno 0 -maf 0.01). After filtering, approximately 220K SNPs were used to phase the genotypes using the software SHAPEIT v.2.r778. 109 For each chromosome separately we ran Chromopainter first to estimate the ''Ne'' and ''mu'' parameters using 10 expectation maximization iteration (-i 10). These parameters were then used to paint each individual against all the others (-a 0 0). Finally we used ''Chromocombine'' to merge the painting information from each chromosome and obtain the normalization parameter ''c''.
The estimated matrix of chunk counts (Data S1C) obtained from Chromocombine was then used as input to the fineSTRUCTURE algorithm. This program was run using 1,000,000 burnin and sampling iterations with sampling every 1000 iterations for the MCMC. Following the method described in Leslie et al. 75 we then extracted the state with the highest posterior probability and performed an additional 100,000 burn-in iterations using the maximum concordance method to obtain the final tree. The information about the optimal number of groups and the cluster assignment of each sample was taken from the file ''.tree'' generated by the program.

Estimated effective migration surface
To visualize how geographical barriers affected migration between populations we used the software EEMS. 66 The same set of nonrelated ancient samples used for the IBD analyses were used to generate a pairwise dissimilarity matrix using the bed2diffs v.2. program. EEMS was initially run using 500 demes with MCMC chains parameters of 100,000 burn-in and 200,000 sampling iterations. This run was repeated 10 times using different random seeds. The run with the highest likelihood was then selected for further refinement using the same number of demes and MCMC settings of 1000,000 burn-in and 2000,000 sampling iterations.

BMI analysis
To investigate the distribution of body mass index across European Neolithic populations we calculated the polygenic risk score (PRS) for 247 individuals using the summary statistics calculated by the Neale Lab (http://www.nealelab.is/uk-biobank) using the UK BioBank resource. Prior to obtaining the PRS information we filtered individuals with more than 30% of BMI SNPs missing. We did not allow missing genotypes to be present in this analysis. SNPs in this dataset were filtered using a clumping/threshold approach through the software PLINK 1.9 with parameters (-clump-p1 0.01 -clump-kb 1000 -clump-r2 0.1). After filtering we obtained approximately 12 thousands SNPs that we used to compute the PRS in 247 ancient samples using the -score option in plink ( Figure S7).  Neolithic, CA=Copper Age, WHG=Western hunter-gatherer, CHG=Caucasus hunter-gatherer, EHG=Eastern hunter-gatherer). The Anatolian Neolithic population is represented by individuals from Marmara, Barcın.

Figure S3: Temporal distribution of hunter-gatherer ancestry in Neolithic
Europeans. Related to Star Methods and Data S2E. qpAdm was used to measure the quantity of Western Hunter-gatherer (WHG) populations, represented by the individuals KO1 and Loschbour, present in European Neolithic populations. Each point represents with its error bar a group with at least 2 individuals from an archaeological site and time period. The WHG percentages are reported with 95% error bars. Figure S4: Comparisons of IBD and F ROH estimates from imputed WGS and SNP captured data. Related to Star Methods and Data S1B, S1D. (A) F ROH were compared between imputations of WGS and SNP captured data (respectively "F ROH WGS" and "F ROH SNPCap") where these are available from the same samples. These correlate with a P-value lower than and the regression line (coloured blue) error margins overlap with the 1:1 plot (red line). Also two yellow coloured points denote where three genomes (NE1, Stuttgart_LBK) also had WGS F ROH estimates available from high coverage SNP calls -these are plotted Vs F ROH SNPCap. (B) Plots of IBD sharing values involving each of nine Neolithic samples for which WGS and SNP capture -based estimates are available. For each, total estimates of the genome shared with all other Neolithic samples is plotted alternately using the two different data sets. The WGS and SNP-derived values correlate significantly (each at p < 10 -15 ) and vary around the 1:1 plotline, drawn in red.  Genotypes imputed from a SNP captured genome (LBK) compared against those called from its high coverage genome sequence using the software SNPSift S43 . Heterozygous and transversion SNPs show overall lower accuracy compared to homozygous and transition imputed SNPs. The SNPs corresponding to the 1240K positions were predicted overall with higher accuracy compared to other variants. (B) Genotypes imputed from a downsampled WGS genome (LBK) also here compared against its high coverage genome sequence. Only transversions are considered in this analysis and heterozygous calls show overall less accuracy than homozygous ones. Overall the percentage of safely imputed genotypes remains high for each genotype class across different probability thresholds.