The ancient Mediterranean port city of Ashkelon, identified as “Philistine” during the Iron Age, underwent a marked cultural change between the Late Bronze and the early Iron Age. It has been long debated whether this change was driven by a substantial movement of people, possibly linked to a larger migration of the so-called “Sea Peoples.” Here, we report genome-wide data of 10 Bronze and Iron Age individuals from Ashkelon. We find that the early Iron Age population was genetically distinct due to a European-related admixture. This genetic signal is no longer detectible in the later Iron Age population. Our results support that a migration event occurred during the Bronze to Iron Age transition in Ashkelon but did not leave a long-lasting genetic signature.

Here, we report genome-wide data from human remains excavated at the ancient seaport of Ashkelon, forming a genetic time series encompassing the Bronze to Iron Age transition ( Fig. 1, A and B ). We find that all three Ashkelon populations derive most of their ancestry from the local Levantine gene pool. The early Iron Age population was distinct in its high genetic affinity to European-derived populations and in the high variation of that affinity, suggesting that a gene flow from a European-related gene pool entered Ashkelon either at the end of the Bronze Age or at the beginning of the Iron Age. Of the available contemporaneous populations, we model the southern European gene pool as the best proxy for this incoming gene flow. Last, we observe that the excess European affinity of the early Iron Age individuals does not persist in the later Iron Age population, suggesting that it had a limited genetic impact on the long-term population structure of the people in Ashkelon.

Recent ancient DNA (aDNA) studies have reported a high degree of genetic continuity in the Levant during the late Pleistocene and early Holocene that was followed by increasing population admixtures with Anatolian- and Iranian-related populations at least up to the Middle Bronze Age ( 11 – 14 ). Genome-wide data from Late Bronze and Iron Age populations have, so far, not been available for this region.

Within the history of the Eastern Mediterranean, the end of the Bronze Age and the beginning of the Iron Age were marked by exceptional cultural disarray that followed the demise of prosperous economies and cultures in Greece, Egypt, the Levant, and Anatolia ( 1 ). During the 12th century BCE, coincident with these events, substantive cultural changes appeared in the archeological record of Ashkelon, Ashdod, and Ekron, three of the five core cities mentioned as “Philistine” in the Hebrew bible ( 2 – 4 ). These settlements were distinct from neighboring sites in architectural traditions and material culture ( 2 – 4 ). Resemblances between the new cultural traits and 13th century patterns found in the Aegean have led some scholars to explain this so-called “Philistine phenomenon” by a migration from an Aegean-related source, potentially associated with the “Sea Peoples,” a population that is thought to have settled in other parts of the coastal Eastern Mediterranean ( 2 , 3 ). This hypothesis has been challenged by those arguing that this cultural change was driven by a diffusion of knowledge or internal development of ideas ( 5 – 8 ) rather than by a large-scale movement of people. Even for those who do accept the idea of large-scale mobility, the homeland of the new arrivals is disputed with suggested alternatives including Cyprus or Cilicia ( 4 ), a mixture of non-Aegean east Mediterranean peoples ( 8 ), and mixed heterogeneous maritime groups, akin to pirates ( 9 ). Proposed links go as far as northern Italy where depopulation events have been suggested to trigger population movements throughout the Mediterranean ( 10 ).

RESULTS

Genome-wide data across the Bronze/Iron Age transition in Ashkelon We extracted and sequenced DNA from 108 skeletal elements excavated in Ashkelon. In line with the low DNA preservation previously reported for the southern Levant (11–14), only 10 yielded sufficient amounts of human DNA (data file S1). Sequencing libraries for these 10 individuals were enriched for human DNA using an established in-solution DNA capture targeting 1.24 million genome-wide single-nucleotide polymorphisms (SNPs) (“1240 k capture”) (15–17). The earliest three individuals were excavated from a Bronze Age necropolis (18–20) and dated by the archeological context to the Middle Bronze IIC–Late Bronze Age II (labeled ASH_LBA; two of the individuals were directly carbon-dated to 1746 to 1542 cal BCE). Four early Iron Age infants were recovered from burials beneath late 12th century Iron Age I houses (labeled ASH_IA1; directly carbon-dated to 1379 to 1126 cal BCE), and three later Iron Age individuals were recovered from a cemetery adjacent to the city wall of ancient Ashkelon, which was estimated to have been used between the 10th and the 9th century BCE (labeled ASH_IA2; one individual directly carbon-dated to 1257 to 1042 cal BCE) (Fig. 1, A and B, Table 1, table S1, and text S1). While the chronological ranges estimated by the archeological context are approximately within the range of the radiocarbon assays conducted directly on the petrous bones sampled for DNA (table S1), we caution that they should be taken as a rough estimate, since the presence of marine carbon in the diet could make the calibrated radiocarbon dates older than the true age. Table 1 An overview of the Ashkelon genomes reported in this study. For each individual, the analysis group is given (ASH_LBA, Ashkelon Late Bronze Age; ASH_IA1, Ashkelon Iron Age 1; ASH_IA2, Ashkelon Iron Age 2). 14C dating results are given in cal BCE in two-sigma range (NA, not available). Detailed dating information is provided in text S1 and table S1. The proportion of human DNA and the mean coverage on 1240 K target sites in the “1240 K” enriched libraries are given. The assigned genetic sex is listed (F, Female; M, Male). Uniparental haplogroups (mt, mitochondrial; Ychr, Y chromosome) are listed. View this table: To control for the quality of the dataset, we estimated exogenous DNA levels and relatedness. Mitochondrial contamination for all 10 individuals ranged between 0 and 9%, with only two of them above 5% (table S2). For males (n = 4), nuclear contamination was estimated to a maximum of 4.3% (table S2). None of the individuals were either first- or second-degree related to each other (fig. S1). We performed population genetic analysis on a merged dataset, including the genotype data of the newly reported individuals from this study and previously published datasets from 638 ancient individuals and 4943 individuals belonging to 298 present-day populations. (11, 12, 21). (An overview of the key ancient individuals analyzed is given in data file S2).

Persistence of the local gene pool in the Bronze Age Levant To understand the genetic profile of the ancient Ashkelon individuals, we began by projecting them onto the first two axes of variation (PC1 and PC2) of present-day west Eurasians, inferred from principal component analysis (PCA) (Fig. 2A and fig. S2) (12). The ASH_LBA individuals overlap with the cline of present-day Near Easterners and are close to earlier Bronze Age Levantines and Anatolians [Early Bronze Age individuals from ‘Ain-Ghazal, Jordan labeled “Jordan_EBA” (12); Middle Bronze Age individuals from Sidon, Lebanon labeled “Lebanon_MBA” (13); and Early to Late Bronze Age individuals from Central Anatolia labeled “Anatolia_EBA” (21) and “Anatolia_MLBA” (21), respectively]. Compared to earlier Levantines [Levantine early farmers from present-day Jordan and Israel labeled “Levant_N” (12) and Chalcolithic individuals from Peqi’in, Israel labeled “Levant_ChL” (11)], the Bronze Age individuals including ASH_LBA are all shifted along PC2 toward ancient Iranian and Caucasus individuals [e.g., the Caucasus Mesolithic hunter-gatherers labeled “CHG” (22); early farmers excavated from present-day Iran labeled “Iran_N” (12); and the Chalcolithic individuals excavated from present-day western Iran labeled “Iran_ChL” (12)], in agreement with previous observations (11–13). Unsupervised genetic clustering using ADMIXTURE (23) shows a similar pattern [with K = 9 clusters; Fig. 2B and fig. S3]; ASH_LBA is assigned the same major ancestral component as all earlier Levantine populations (shown in orange). Consistent with their PCA positions, a second component (shown in green) that is maximized in Iran_N is, on average, higher in ASH_LBA and in each of the earlier Bronze Age Levantines compared to all earlier Levantines (20 to 30% and 3 to 8%, respectively). Fig. 2 PCA and ADMIXTURE analysis. (A) Ancient genomes (marked with color-filled symbols) projected onto the principal components inferred from present-day west Eurasians (gray circles) (fig. S2). The newly reported Ashkelon populations are annotated in the upper corner. (B) ADMIXTURE analysis. A selected set of ancient individuals (as well as present-day Sardinians) is plotted (K = 9 was chosen since it is the cluster number that maximizes components correlated to the most differentiated populations in the west Eurasian PCA). To formally test the qualitative observations based on the PCA and ADMIXTURE analyses, we used f4-statistics (24). Consistent with the ASH_LBA positions on PC2, the f4-statistic of the form f4 (ASH_LBA, Levant_N/Levant_ChL; test, Mbuti) (figs. S4 and S5) estimated excess allele sharing between ASH_LBA and ancient Iranian/Caucasus-related populations (such as CHG, Iran_N, and Iran_ChL) compared to the Neolithic/Chalcolithic Levantines (Levant_N and Levant_ChL), confirming previous reports of post-Neolithic gene flows from prehistoric populations related to Iran or the Caucasus into the Levant (11–13). Accordingly, modeling ASH_LBA as a two-way admixture using qpAdm (12, 17) produced a fitting model (χ2P = 0.445), in which ASH_LBA derives around 60% of their ancestry from Levant_ChL (60.0 ± 6.5%; estimate ±1 SE) and the rest from Iran_ChL. The spatiotemporal origins of this post-Neolithic gene flow remain broadly defined since alternative two-way models also fit, albeit with smaller P values. In these models, Anatolia_EBA/Anatolia_BA/Armenia_ChL are used as the second source, replacing Iran_ChL (χ2P = 0.053 to 0.060; table S3). We then examined whether there are noticeable genetic differences between ASH_LBA and the earlier Bronze Age Levant populations by measuring f4 (ASH_LBA, Jordan_EBA/Lebanon_MBA; test, Mbuti (figs. S6 and S7). We find that ASH_LBA and both earlier Bronze Age populations are symmetrically related to all test populations within our data’s resolution (|f 4 | < 3 SE). However, we observe a marginal excess of affinity between ASH_LBA and populations genetically related to ancient Iran and the Caucasus (such as CHG and Steppe_Eneolithic) when compared to Jordan_EBA (fig. S7). Within our current resolution, we cannot distinguish whether this affinity is due to a small-scale gene flow from an Iranian/Caucasus-related population entering the Levant between the Early and Late Bronze Age or whether it reflects a certain population structure during this period. Nonetheless, the apparent symmetry suggests a high degree of genetic continuity throughout at least a millennium between culturally and geographically distinct Bronze Age Levantine groups occupying a region that stretches from the inland southern Levant where present-day Jordan is located and along the coastal regions of present-day Israel and Lebanon.

Genetic discontinuity between the Bronze Age and the early Iron Age people of Ashkelon In comparison to ASH_LBA, the four ASH_IA1 individuals from the following Iron Age I period are, on average, shifted along PC1 toward the European cline and are more spread out along PC1, overlapping with ASH_LBA on one extreme and with the Greek Late Bronze Age “S_Greece_LBA” (25) on the other (Fig. 2A). Similarly, genetic clustering assigns ASH_IA1 with an average of 14% contribution from a cluster maximized in the Mesolithic European hunter-gatherers labeled “WHG” (shown in blue in Fig. 2B) (15, 22, 26). This component is inferred only in small proportions in earlier Bronze Age Levantine populations (2 to 9%). In light of these observations, we formally tested the variation within each of the three Ashkelon populations by computing the variance of the square rooted statistic f4 (Ashkelon individual 1, Ashkelon individual 2; test, Mbuti) (fig. S8). Consistent with the wide spread of ASH_IA1 in the west Eurasian PCA, the variance is significantly higher in ASH_IA1 than in either ASH_LBA or ASH_IA2 (P < 2.2 × 10−16 for both by a Fligner-Killeen test) (table S4 and fig. S8), suggesting that the ASH_IA1 individuals were more heterogeneous in their genetic affinities to global populations in respect to both ASH_LBA and ASH_IA2. We next formally quantified the genetic differences between ASH_IA1 and the earlier ASH_LBA by calculating f4 (ASH_IA1, ASH_LBA; test, Mbuti) (Fig. 3A). In agreement with the PCA and ADMIXTURE results, only European hunter-gatherers (including WHG) and populations sharing a history of genetic admixture with European hunter-gatherers (e.g., as European Neolithic and post-Neolithic populations) produced significantly positive f4-statistics (Z ≥ 3), suggesting that, compared to ASH_LBA, ASH_IA1 has additional European-related ancestry. To estimate the levels of the European-related ancestry in all Ashkelon populations, we compared qpAdm models of two-way mixtures (Levant_ChL and Iran_ChL; i.e., 0% contribution from WHG) to three-way ones, in which we add WHG as the third source (Fig. 3C and table S5). ASH_LBA and ASH_IA2 fit well with the two-way model (χ2P = 0.445 and χ2P = 0.313, respectively; table S5), whereas the three-way one infers small nonsignificant proportions (−2.3 ± 2.2 and 2.0 ± 2.2% for ASH_LBA and ASH_IA2, respectively; table S5). In contrast, for ASH_IA1, the two-way model is inadequate (χ2P = 3.80 × 106; table S5) and the three-way one fits (χ2P = 0.765). Thus, of the three, additional WHG-like ancestry is only necessary to model ASH_IA1. Fig. 3 European-related admixture detected in ASH_IA1. (A) ASH_IA1 shares access affinity with European-related populations compared to ASH_LBA. We plot the top and bottom 40 values of f4 (ASH_IA1, ASH_IA2; X, Mbuti) on the map. Circles mark the ancient populations and triangles the present-day ones. Z-scores calculated by 5-centimorgan block jackknifing are represented by the size of the symbols. “X” share more alleles with ASH_IA1 when values are positive and with ASH_IA2 when negative. The five groups with the most positive values are annotated on the map (Z > 2.3). (B) We plot the ancestral proportions of the Ashkelon individuals inferred by qpAdm using Iran_ChL, Levant_ChL, and WHG as sources ±1 SEs. P values are annotated under each model. In cases when the three-way model failed (χ2P < 0.05), we plot the fitting two-way model. The WHG ancestry is necessary only in ASH_IA1. We find that the PC1 coordinates positively correlate with the proportion of WHG ancestry modeled in the Ashkelon individuals (fig. S9 and table S6), suggesting that WHG reasonably tag a European-related ancestral component within the ASH_IA1 individuals. However, these Mesolithic individuals are unlikely to be a good proxy for the true source in the much later early Iron Age. To examine more proximate sources, we compiled a set of chronologically and geographically relevant candidate populations, including populations that shared higher affinity with ASH_IA1 compared to ASH_LBA in the above f4-statistic. Subsequently, we modeled ASH_IA1 as two-way and three-way mixtures of ASH_LBA and combinations of the candidate populations (table S7). Of the 51 tested models, we find four plausible ones (χ2P > 0.05), all are two-way mixtures. The best supported one (χ2P = 0.675) infers that ASH_IA1 derives around 43% of ancestry from the Greek Bronze Age “Crete_Odigitria_BA” (43.1 ± 19.2%) and the rest from the ASH_LBA population. ASH_IA1 could also be modeled with either the modern “Sardinian” (35.2 ± 17.4%; χ2P = 0.070), the Bronze Age “Iberia_BA” (21.8 ± 21.1%; χ2P = 0.205), or the Bronze Age “Steppe_MLBA” (15.7 ± 9.1%; χ2P = 0.050) as the second source population to ASH_LBA. To check whether these results are due to the low coverage of ASH_LBA, we repeated this analysis, but this time, we modeled ASH_IA1 as a three-way mixture of each of the candidate populations, Levant_ChL and Iran_ChL. The two latter populations have higher genome coverage and can model ASH_LBA well in combination (table S3). In this analysis, only the models including “Sardinian,” “Crete_Odigitria_BA,” or “Iberia_BA” as the candidate population provided a good fit (χ2P = 0.715, 49.3 ± 8.5%; χ2P = 0.972, 38.0 ± 22.0%; and χ2P = 0.964, 25.8 ± 9.3%, respectively). We note that, because of geographical and temporal sampling gaps, populations that potentially contributed the “European-related” admixture in ASH_IA1 could be missing from the dataset. Therefore, better proxies might be found in the future when more data is available. Nonetheless, the tested candidate populations from Anatolia, Egypt, and the Levant that did not produce well-fitting models can be excluded as potential sources of the admixture observed in ASH_IA1.