Ancient genomics reveals tripartite origins of Japanese populations

Our initial screening focused on 14 ancient skeletal remains excavated from six archaeological sites across the archipelago (see note S1). High levels of endogenous human DNA were preserved in 12 of these samples (Table 1), which were then further shotgun-sequenced to higher coverage, ranging from 0.88× to 7.51× (Fig. 1 and table S1). Nine of the 12 samples are associated with the Jomon culture, representing the west and central parts of the archipelago and four different stages of the Jomon period (Initial, Early, Middle, and Late Jomon) (Fig. 1). The remaining three samples date to ~1.3 ka ago, falling within the Kofun period. We confirm that all newly sequenced genomes show postmortem damage patterns (fig. S1) and a low level of modern human contamination (<2.15%) (Table 1 and table S2). Our kinship analysis confirms that all pairs of individuals are unrelated (fig. S2). Mitochondrial haplogroups for all Jomon individuals belong to the N9b or M7a clades, which are strongly associated with this population (11–14, 22) and rare outside of Japan today (23). The three Jomon males (table S3) belong to the Y chromosome haplogroup D1b1, which is present in modern Japanese populations but almost absent in other East Asians (24). In contrast, the Kofun individuals all belong to mitochondrial haplogroups that are common in present-day East Asians (25), while the single Kofun male has the O3a2c Y chromosome haplogroup, which is also found throughout East Asia, particularly in mainland China (26). To place our data within the wider context of East Eurasian demography, we combined the ancient Japanese genomes with genomic data from previously published ancient (table S4 and fig. S3) and present-day individuals. Throughout this study, the modern Japanese population is represented by either data from the Simons Genome Diversity Project (SGDP) (27) or JPT (i.e., Japanese in Tokyo) in the 1000 Genome Project phase 3 (28). However, we note that ancestral heterogeneity exists across the archipelago today, which is not fully captured by this standard reference set. Other ancient and present-day populations analyzed in this study are primarily labeled by either geographic or cultural contexts.

We explored the genetic diversity within our time series data by looking at the shared genetic drift between all pairwise comparisons of individuals from both the ancient and modern (SGDP) Japanese populations using the statistic f₃(Individual_1, Individual_2; Mbuti) (Fig. 2A). Our results very clearly define three distinct clusters of Jomon, Yayoi, and Kofun individuals, the last of whom group with the modern Japanese individuals, suggesting that cultural shifts were accompanied by genomic changes. Despite the large spatial and temporal variation in the Jomon dataset, extremely high levels of shared drift are observed between all 12 individuals. The Yayoi individuals are most closely related to each other and also have a higher affinity to the Jomon than to the Kofun individuals. The Kofun and modern Japanese individuals are almost indistinguishable from one another by this metric, implying some level of genetic continuity over the past 1400 years.

We further investigated the genome-wide autosomal affinities of our ancient Japanese individuals to continental populations using principal components analysis (PCA). We projected ancient individuals onto the genetic variation of present-day populations in the SGDP dataset from South and Central Asia, Southeast and East Asia, and Siberia (Fig. 2B and fig. S4). We observe that the ancient Japanese individuals separate into their respective cultural designations along PC1. All Jomon individuals form a tight cluster, falling apart from other ancient populations, as well as present-day Southeast and East Asians, suggesting a sustained geographic isolation. The two Yayoi individuals appear near to this Jomon cluster, supporting genetic and morphological similarity to Jomon as reported in (15, 16). However, a shift toward East Asian populations implies the presence of additional continental ancestry in Yayoi. Ancient individuals from East Eurasia show a geographic cline from the south to the north on PC2: Southern China, Yellow River, Northern China, West Liao River, Devil’s Gate Cave, Amur River, and Baikal. The three individuals from the Kofun period fall within the diversity of the Yellow River cluster.

ADMIXTURE analysis with the Human Origins Array dataset also supports the increase in continental gene flow into the archipelago after the end of the Jomon period (Fig. 2C and figs. S5 and S6). The Jomon have a distinct ancestral component (represented in red in Fig. 2C), which is also visible at high levels in the Yayoi and remains at reduced levels in the Kofun and Japanese. New ancestral components appear in the Yayoi, at proportions similar to the profile seen in the Amur River basin and surrounding regions. These include a larger component that is dominant in Northeast Asians (represented by light blue) and another smaller component that represents much broader East Asian ancestry (represented by yellow). This East Asian component becomes dominant in the Kofun period and the modern Japanese population.

The separation of the Jomon from other populations (Fig. 2) supports the idea that they form a distinct lineage among East Eurasians, as proposed in previous studies (13, 14). To explore the depth of this divergence, we reconstructed the phylogenetic relationship of Jomon with 17 other ancient and present-day populations using TreeMix with differing numbers of admixture events (Fig. 3A and fig. S7) (29). Our results infer that Jomon emerged after the early divergences of Upper Paleolithic East Eurasians (Tianyuan and Salkhit) and ancient Southeast Asian hunter-gatherers (Hoabinhian), but before the splitting off of other samples including present-day East Asians, an ancient Nepali (Chokhopani), hunter-gatherers from Baikal (Shamanka_EN and Lokomotiv_EN) and Chertovy Vorota Cave (Devil’s Gate Cave) in the Primorye Region, and a Pleistocene Alaskan (USR1). We further confirm the position of Jomon between two other deeply diverged hunter-gatherer lineages in this tree through formal tests of symmetric models using f₄(Mbuti, X; Hoabinhian/DevilsCave_N, Jomon) (fig. S8, A and B). These show that all East Asian individuals in our dataset since the beginning of the Jomon period have a higher affinity to Jomon than the earlier diverged Hoabinhian but have a lower affinity in comparison to DevilsCave_N; this supports the inferred phylogeny of three distinct hunter-gatherer lineages in East Eurasia, rather than a previously proposed model in which Jomon is a mixture of Hoabinhian and East Asian-related lineages [f₃(Jomon; Hoabinhian, DevilsCave_N) = 0.193, Z = 61.355; table S5] (12, 30). We also consistently infer gene flow from Jomon to the modern Japanese across all migration models tested, with genetic contributions ranging from 8.9 to 11.5% (fig. S7). This is consistent with the mean Jomon component of 9.31% in the present-day Japanese individuals estimated from our ADMIXTURE analysis (Fig. 2C). These results suggest a deep divergence of the Jomon and an ancestral link to present-day Japanese populations.

We apply population genetic modeling to estimate the timing of the appearance of the Jomon lineage. Our approach makes use of genome-wide patterns of runs of homozygosity (ROHs) to identify the demographic scenario that best fits the ROH spectrum observed in our oldest and highest coverage sample, JpKa6904 (see note S2). The distribution of ROH tracts reflects the effective population size and the time to the most recent common ancestor between two copies of haplotypes within an individual (31, 32). The 8.8-ka-old Jomon carries high levels of ROH, in particular with the highest frequency of short ROH (due to population effects rather than recent inbreeding) yet reported (Fig. 3B) (33). This pattern, coupled with strong shared genetic drift among the Jomon individuals (fig. S9), implies that the Jomon population underwent a severe population bottleneck. With a search of parameter space of population size and split time, our estimates place the appearance of the Jomon lineage between 15 and 20 ka ago, followed by the maintenance of a very small population size of ~1000 until at least the Initial Jomon period (Fig. 3C and figs. S10 and S11). This coincides with rising sea levels and the severing of the land bridge to the mainland at the end of the LGM (34, 35) and shortly precedes the first appearance of Jomon pottery in the archipelago (2).

We then asked whether the Jomon had any contact with continental Upper Paleolithic people after the divergence of their lineage, but before their isolation in the archipelago, using the statistic f₄(Mbuti, X; Jomon, Han/Dai/Japanese) (fig. S8, C to E). Among the Upper Paleolithic individuals tested, only Yana_UP is significantly closer to Jomon than Han, Dai, or Japanese, respectively (Z > 3.366). This affinity is still detectable even if we replace these reference populations with the other Southeast and East Asians (table S6), supporting gene flow between the ancestors of Jomon and Ancient North Siberians, a population widespread in North Eurasia before the LGM (19).

We lastly investigated potential temporal and spatial variation within the Jomon population. Three temporal groups, defined by the Initial, Early, and Middle-Late-Final stages of the Jomon period, show similar levels of shared genetic drift with ancient and present-day continental populations, implying little or no genetic impact from outside of the archipelago over these subperiods (Fig. 3D and fig. S12). This pattern is further backed up by the absence of any significant gene flow, observed with the statistic f₄(Mbuti, X; sub_Jomon_i, sub_Jomon_j), where i and j are any pairs of the three Jomon groups (fig. S13). These Jomon individuals similarly do not display any diversity in genetic affinity to continental populations when grouped by geography (i.e., different islands on which the samples are located: Honshu, Shikoku, and Rebun Island) (fig. S14). The only observable difference within the Jomon individuals is a slightly higher affinity between sites located on Honshu, implying an insular effect with restricted gene flow between Honshu and other islands (fig. S15). Overall, these results show limited spatial and temporal genetic variation within the Jomon population, supporting the idea of near-complete isolation from the rest of Asia for thousands of years.

Two individuals from the southwestern part of the archipelago associated with the Yayoi culture (Fig. 1) (15) were found to have both Jomon and continental ancestry (Fig. 2). Our qpAdm analysis rejected a model in which the Yayoi are unadmixed descendants of the Jomon (P = 0.00003), as opposed to a morphological assessment that groups these two individuals as a part of the Jomon lineage (16). The non-Jomon ancestral component could have been introduced by people who brought rice cultivation to the archipelago. We first tested whether any ancient East Eurasian population has higher genetic affinity to the Yayoi than Jomon using f₄(Mbuti, X; Jomon, Yayoi) (Fig. 4A and fig. S16). Most of the sampled ancient populations in the continent do not show significant affinity to the Yayoi, including populations from the Yellow River basin (20) where rice farming first spread from the lower Yangtze Valley (a hypothetical origin of japonica rice, i.e., wet rice) (36). However, excess affinity to the Yayoi was detected in populations who had no cultural relation to rice farming (Z > 3.0): those from the West Liao River basin in Northeast China (WLR_BA_o and HMMH_MN), Baikal (Lokomotive_EN, Shamanka_EN, and UstBelaya_EBA), and Northeast Siberia (Ekven_IA). This affinity is not observable in two more individuals (WLR_BA; Z = 1.493) who came from the same archaeological site as the other Bronze Age West Liao River individual (WLR_BA_o) (20). These two individuals have much higher Yellow River–related ancestry (81.4 ± 6.7%) than WLR_BA_o (1.8 ± 9.1%) (20), which implies that the ancient Yellow River populations tested are unlikely to be a major source of the non-Jomon ancestry carried by the Yayoi.

To further distinguish these six potential sources of continental ancestry, we modeled the Yayoi as a two-way admixture of the Jomon and each in turn using qpWave (table S7). The admixture model was confidently supported (P > 0.05) for three of these: Baikal hunter-gatherers and the West Liao Middle Neolithic or Bronze Age individuals with a high level of Amur River ancestry (20). These groups all share a dominant Northeast Asian ancestral component (Fig. 2C and fig. S17) (20). Admixture fractions of 55.0 ± 10.1%, 50.6 ± 8.8%, or 58.4 ± 7.6% for the Jomon input were estimated by qpAdm when each of these three were, respectively, used as the second source (Fig. 4B and table S7), with a fraction of 61.3 ± 7.4% returned when the West Liao Middle Neolithic and Bronze Age individuals were merged into a single source population (table S7).

We further confirm by f₄(Mbuti, Jomon; Yayoi_1, Yayoi_2) that the level of Jomon ancestry is comparable between these two Yayoi individuals (Z = 1.309). These results imply an approximately equal ratio of indigenous hunter-gatherer and migrant contribution to the Yayoi communities associated with this northwestern Kyushu site. This parity is particularly notable when compared to agricultural migrations in western Eurasia, where minimal hunter-gatherer contribution is observed in many regions, including the archipelago of Britain and Ireland (37–40), which mirrors Japan as an insular geographic extreme of Eurasia.

While the West Liao populations used in our admixture models did not themselves practice rice farming, they are situated just north of a hypothesized route of agricultural spread into Japan, to which our results lend weight. This follows the Shandong Peninsula (northeastern China) into the Liaodong Peninsula (northwestern part of the Korean Peninsula) and then reaches the archipelago via the Korean Peninsula (41).

We further investigated how the Yayoi culture spread into the archipelago using outgroup f₃ statistics that measured genetic affinity between the Yayoi and each of the Jomon individuals. We find that the strength of shared genetic drift has a significant correlation with the distance from the location of the Yayoi individuals (P = 0.00697; Fig. 4C); the closer the Jomon archaeological sites are to the Yayoi site, the more the Jomon individuals share genetic drift with the Yayoi. This result supports the introduction of rice farming via the Korean Peninsula (41, 42), followed by admixture with local Jomon populations in the south of the archipelago.

Historical records provide strong support for continued population movement from the continent to the archipelago during the Kofun period (1). However, our qpWave modeling of the three Kofun individuals rejected the two-way admixture of Jomon and Northeast Asian ancestry that fitted the Yayoi individuals (P < 0.05; table S8). Thus, the Kofun are genetically distinct from the Yayoi in terms of their ancestral components, as supported from our outgroup f₃, PCA, and ADMIXTURE clusters (Fig. 2), as well as from previous morphological studies (43, 44). To identify additional ancestral groups that contributed to the genetic makeup of the Kofun individuals, we tested genetic affinity between the Kofun and each of the continental populations using f₄(Mbuti, X; Yayoi, Kofun) (Fig. 5A and fig. S18). We find that most of the ancient or modern populations in our dataset are significantly closer to the Kofun than they are to the Yayoi. This finding implies additional migration to the archipelago during the six centuries that separates the genomes from these two periods.

We attempted to narrow the source of this migration by testing the fit of two-way admixtures of the Yayoi and populations identified from our f₄ statistics as being significantly closer to the Kofun (table S9). This admixture model was confidently supported with P > 0.05 for only 5 of 59 populations tested. We then applied qpAdm to quantify the genetic contribution from Yayoi and each of these sources in turn (table S9). The two-way admixture model was subsequently rejected in an additional two populations due to a lack of support from different reference sets. The remaining three populations (Han, Korean, and YR_LBIA) show 20 to 30% contributions to the Kofun (fig. S19). These three all share strong genetic drift (highlighted by the green square in fig. S17), characterized with a major component of broadly East Asian ancestry in their ADMIXTURE profiles (fig. S6). To further screen the source of additional ancestry in the Kofun individuals, we tested a three-way admixture by replacing the Yayoi ancestry with Jomon and Northeast Asian ancestry (table S10). Only Han were successfully modeled as a source of ancestry in the model (Fig. 5B), with a significantly better fitting of the three-way admixture than any possible two-way admixture models (table S11). Given that Jomon ancestry is diluted by approximately a factor of 4 between the Yayoi and Kofun populations sampled, these results suggest that the state formation phase saw a large influx of migrants who had East Asian ancestry.

We then explore the possibility that the continental ancestry observed in both the Yayoi and Kofun periods derives from the same source, with intermediate levels of Northeast and East Asian ancestry (table S12). Only one candidate was found to better fit a two-way mixture for Kofun, a population of the Late Bronze Age and Iron Age individuals from the Yellow River basin (YR_LBIA) (20), although this was not consistent across the reference sets (nested, P = 0.100; table S13). Despite not showing statistically significant gene flow with Yayoi to the exclusion of Jomon (Z = 2.487) (Fig. 4A and fig. S16), we find that a two-way model between YR_LBIA and Jomon also fits the Yayoi (table S14). This Yellow River population has an intermediate genetic profile with approximately 40% Northeast Asian and 60% East Asian (i.e., Han) ancestry, as estimated by qpAdm. Thus, this is an intermediate genetic profile that fits both Yayoi and Kofun in certain models, with an input of 37.4 ± 1.9% and 87.5 ± 0.8% in these populations (table S14). These results imply that continuous gene flow from a single source may be sufficient to explain the genetic changes between the Yayoi and Kofun.

However, our broader analyses strongly suggest that a single source of gene flow is less likely than two distinct waves of migrations. First, the ratios of Northeast Asian to East Asian ancestry identified in ADMIXTURE were starkly different between the Yayoi (1.9:1) and the Kofun (1:2.5) (Fig. 2C). Second, this contrast in continental affinity is also observable in the different forms of f statistics, in which there is a repeated pattern of Yayoi having significant affinity to those with Northeast Asian ancestry (Fig. 4A and fig. S16), while the Kofun individuals form a tight cluster with other East Asians including Han and ancient Yellow River populations (Fig. 5A and figs. S17 and S18). Last, we find support for a two-pulse model from our dating of the admixture in the Kofun individuals by DATES (fig. S20) (45). A single admixture event with the intermediate population (i.e., YR_LBIA) is estimated to have occurred 1840 ± 213 years before the present (B.P.), which is much later than the onset of the Yayoi period (~3 ka ago). In contrast, if two separate admixture events with two distinct sources are assumed, the resulting estimates reasonably fit the timings consistent with the beginning of the Yayoi and Kofun periods (3448 ± 825 years B.P. for the admixture between Jomon and Northeast Asian ancestry and 1748 ± 175 years B.P. for Jomon and East Asian ancestry; fig. S20). These genetic findings are further supported by both the archaeological evidence and the historical records, which document the arrival of new people from the continent during the period (1).

The three Kofun individuals are genetically similar to the present-day Japanese, as shown in Fig. 2. This implies no substantial change in the genetic makeup of Japanese populations since the Kofun period. To look for signals of additional genetic ancestry in the present-day Japanese samples, we tested whether the continental populations have preferential affinity to modern genomes relative to Kofun using f₄(Mbuti, X; Kofun, Japanese) (fig. S21). Although some of the ancient populations exhibit higher affinity to Kofun than Japanese, none of them is supported by qpAdm as an additional source of ancestry present in the Kofun (see note S3 and table S15). Unexpectedly, no ancient or modern population shows additional gene flow with the present-day Japanese to the exclusion of the Kofun. Our admixture modeling further confirms that the present-day Japanese population is sufficiently explained by Kofun ancestry without increasing Jomon or Yayoi ancestry or without introducing an additional ancestor represented by present-day Southeast or East Asians or Siberians (table S16). We also find that the modern Japanese population has the same set of ancestral components as the three-way admixture in the Kofun individuals (table S17), with a slightly increased level of East Asian ancestry in the modern Japanese compared with the Kofun individuals sampled (fig. S22). This suggests a certain level of genetic continuity but not absolute. A strict model of continuity between the Kofun and the modern Japanese population (i.e., with no genetic drift unique to the Kofun lineage) is rejected (table S18) (46). However, we note no dilution of Jomon ancestry in the Japanese population (15.0 ± 3.8%), relative to the Kofun individuals (13.1 ± 3.5%) (fig. S22), as opposed to the preceding Yayoi and Kofun periods where Jomon ancestry significantly decreased due to the continental migrations (Figs. 4 and 5). Testing the genetic cladality between Kofun and Japanese by qpAdm with a “no admixture” model, we found that Kofun forms a clade with Japanese (P = 0.769). These results suggest that the genetic profile of three major ancestral components established by the state formation period has become a foundation for present-day Japanese populations, as also supported from dental and nonmetrical cranial traits (47, 48).

Petrous bones and a tooth from a total of 14 ancient Japanese individuals (11 Jomon and 3 Kofun) from six archaeological sites across the west and central parts of the Japanese archipelago were sampled. Detailed information of each site is provided in note S1, including archaeological contexts and radiocarbon dates that are taken from OxCal4.4 (58) with the IntCal20 curve.

Sample processing was carried out in dedicated ancient DNA facilities at Kanazawa University and Trinity College Dublin where all the precautions for ancient DNA processing were followed as described in (59). Samples were photographed extensively before further processing. All bones were exposed to ultraviolet light for 15 min on either side to remove surface contaminants. Further cleaning of their surface was performed with a drill before extraction. A triangular wedge section of the otic capsule region of the petrous temporal bone was targeted for sampling. Half of the wedge or fragmented tooth was pulverized for DNA extraction. An aliquot of ~0.1 g of the bone powder was subjected to a silica column method with or without an initial washing step by 0.5% bleach solution as in (60) or by predigestion solution as in (61). DNA extracts were purified with MinElute silica columns (QIAGEN) and eluted at a volume of 55 μl of elution buffer.

The initial screening of each sample and blank controls was performed by constructing a double-stranded DNA next-generation sequencing library mainly using the method outlined in (62) with modifications as in (38) but also using a NEBNext Ultra II DNA Library Prep Kit for Illumina with bead-based size selection [<100 base pairs (bp)]. Every library was screened on a MiSeq Illumina platform. DNA extracts with human endogenous content of >20% were selected for high-coverage sequencing and subsequently treated with uracil-DNA-glycosylase (UDG) (63) before a second library construction with the former method. To increase complexity when sequencing for high-coverage, several polymerase chain reactions (PCRs) with unique indexes were prepared from each UDG-treated library. Different PCR cycles were used according to the template concentration. The concentration and quality were then assessed using the Agilent TapeStation 2200 system with the D1000 ScreenTape. High-coverage sequencing for UDG-treated libraries from 12 individuals was carried out on a HiSeq 2500 or NovaSeq 6000 Illumina platform (100-bp single-end reads or 50-bp paired-end reads) at Macrogen (Republic of Korea) or TrinSeq (Ireland).

We trimmed adapters from raw single-end sequence reads using cutadapt v1.9.1 (64) with a minimum overlap of 1 bp and removed sequences shorter than 34 bp in length and from raw paired-end sequence reads using AdapterRemoval v2.2.2 (65) with the following parameters: --collapse --minlength 25 --minadapteroverlap 1 --minquality 25--trimns and --trimqualities. The trimmed reads were aligned to the hg19 reference genome with the revised Cambridge Reference Sequence (rCRS) for mitochondrial DNA (mtDNA) using bwa-aln in BWA (Burrows-Wheeler Aligner) v0.7.5 (66) with relaxed parameters “-l 16500 -n 0.01 -o 2.” We applied a mapping quality of 20 to the aligned data and removed PCR duplicates using SAMtools v1.7 (67). The authenticity of ancient DNA was determined in non–UDG-treated libraries by looking at degradation patterns at 5′- and 3′-end, respectively, using mapDamage2.0 (68). We added read groups using Picard tools version 1.101 (http://broadinstitute.github.io/picard/), before all libraries were merged into one bam per individual. Two final processing steps were carried out on the merged files: Indels were realigned using Genome Analysis Toolkit (GATK) version 3.7-0 (69), and the quality of the two bases at the end of each read was manually reduced to a score of 2 (“softclipped”), to avoid calling genotype from these damage-prone regions. We used the same processing pipeline for published ancient genomes in FASTQ format with the exception of using AdapterRemoval version 2.2 (65) instead of cutadapt; in some cases, BAM files were realigned to our reference genome and then put through the same pipeline. Quality assessment of raw data was conducted with FastQC v0.11.4 (70).

To determine the mitochondrial contamination rate and assign a mitochondrial haplogroup, we aligned sequence reads to the rCRS mitochondrial genome and reprocessed the aligned data using the same pipeline described in the previous section. Consensus sequences were determined to a minimum depth of consensus of 5 and a base quality of 30 using mpileup and bcftools within the SAMtools package. We then assigned a specific haplotype to each consensus sequence using HAPLOFIND (71).

We estimated contamination rates as previously described in (33): Briefly, we calculated the rate of secondary bases that did not match the consensus sequence at haplotype-defining and private mutations defined by HAPLOFIND. We also report the rate when removing single-nucleotide polymorphisms (SNPs) marked as C or G that may be due to postmortem damage (table S2).

To determine the sex of each of our ancient samples, we filtered for reads with a minimum mapping quality of 30 and followed the method outlined in (72). In summary, we calculated the ratio of Y chromosome reads to the total number of sex chromosome (R_y) (72) and assigned female if R_y < 0.016 or male if R_y > 0.075 (table S3). The kinship estimation for all ancient individuals was done using “Relationship Estimation from Ancient DNA” (READ) (73). Any first-degree relatives were removed from analysis.

All newly sequenced individuals determined to be male were piled-up using GATK version 3.7-0 (69) and manually compared to the International Society of Genetic Genealogy (ISOGG) SNP index with a mapping quality of 20 and a base quality of 30.

We created a large panel of ancient East Eurasians that include ancient Japanese newly sequenced in this study and published ancient genomes from Central and Eastern Steppe, Siberia, Southeast Asia, and East Asia (see a full list of ancient individuals in table S4). We genotyped all ancient individuals for the biallelic SNP sites present on two different reference panels in which they were merged to the following: (i) the Human Origin Array (HOA) consisting of 594,896 SNP sites for 1963 modern, ancient, and reference genomes (39); and (ii) the SGDP consisting of 278 modern, ancient, and reference genomes (27) that have been filtered for autosomal transversion-only SNPs with a minor allele frequency of 1%, which left 3,867,366 SNP sites in our merged data. For each SNP site, we randomly called a high-quality single base (bq30) per position to create a pseudo-diploid genotype using GATK version 3.7-0 (69).

PCA was conducted using the smartpca (v16000) from the EIGENSOFT package (v7.2.0) (74) We projected all ancient Japanese and a subsection of other ancient individuals onto present-day East Eurasians in the SGDP panel (n = 112), filtered for transversions and global minor allele frequencies of 1%, with the options of “killr2: YES,” “r2thresh: 0.2,” “numoutlieriter: 0,” “lsqproject: YES,” and “autoshrink: YES.” With the exception of “Yayoi_1,” we included ancient individuals with a coverage of at least 100,000 SNPs.

We used ADMIXTURE v.1.3.0 (75) for unsupervised genetic clustering of ancient East Eurasian samples with at least 100,000 SNPs in the SGDP dataset (n = 189) and present-day populations from the HOA (n = 786; Mbuti, Sardinian, South Asians, Southeast Asians, East Asians, Siberians, Oceanians, and Native Americans). SNPs were pruned for sites in linkage disequilibrium using PLINK v1.90b4.4 (76), with a sliding window size of 50 variants, a step size of 5 variants, and an r² threshold of 0.2 (--indep-pairwise 50 5 0.2), leaving 186,856 SNPs for analysis. We ran 10 replicates with random seeds for the number of clusters (K) from 2 to 12 and chose the run with the minimum cross-validation error for plotting.

Maximum likelihood trees were inferred under various admixture models using TreeMix (v1.13) (29). A subset of ancient and present-day populations in the SGDP dataset was chosen to represent different ancestors across East Eurasia: Ust_Ishim, Yana_UP, MA1, Tianyuan, Salkhit, Papuan, Hoabinhian (La368 and Ma911), Kusunda, Jomon, Chokhopani, Shamanka_EN, Lokomotiv_EN, DevilsCave_N, USR1, Han, Ami, and Japanese. The dataset was filtered for nonmissingness, which left 98,687 SNPs for analysis. We used Mbuti as an outgroup to root the tree and fitted models with no-migration or migrations from one to five to the data. Total 1000 to 1500 replicates were run for each model, of which the tree with a likelihood of the closest to a mean across the replicates was considered the most representative under a given model.

We calculated f statistics using qp3Pop (v300) and qpDstat (v662) with the f₄ mode in the AdmixTools v6.0 package (77). We used outgroup f₃ statistics to measure a genetic relationship between two populations or between all pairs of individuals within a population by specifying Mbuti as an outgroup. SEs in f₄ statistics were calculated from the default block jackknife approach.

We used HaplotypeCaller, part of GATK version 3.7.0, to call diploid genotypes for the oldest Jomon, JpKa6904, and other published high-coverage ancient genomes, including Yana1 (19), USR1 (78), Loschbour, Stuttgart_LBK (39), Mesolithic Irish (SRA62), and Neolithic Irish (JP14) (33). Our diploid call was based on the 1000 Genome phase 3 release v5 panel (28) filtered for transversion sites and a minor allele frequency of 1%. We used VCFtools v0.1.13 (79) to filter our results for a genotype quality of 30 and a minimum depth of 10. Transversion SNPs common to each diploid ancient sample (total of 660,027 SNPs) were used for this analysis. ROHs were measured in each diploid genome using PLINK v1.90b4.4 (76) with the following options: -homozyg --homozyg-density 50 --homozyg-gap 100 --homozyg-kb 500 --homozyg-snp 50 --homozyg-window-het 1 --homozyg-window-snp 50 --homozyg-window-threshold 0.05.

To estimate the population size (N) and split time (T) of the Jomon lineage, we applied coalescent simulations of the Out-of-Africa model (80) to approximating a likelihood surface by fitting genome-wide patterns of ROHs to those observed in our oldest Jomon individual, JpKa6904 (see note S2). Briefly, we first measured the proportion of the Jomon genome under ROH using PLINK v1.90b4.4 (76). Second, we conducted coalescent simulations under different combinations of N and T using ms (81), with an option of “-eA” that allows sampling of ancient chromosomes at a given time in the past. We then fitted the spectrum of ROH fragments ranging from 0.5 to 100 Mb.

We broadly searched for the parameter space from all combinations of 500 ≦ N ≦ 2500 and 10 ≦ T ≦ 40 ka ago using the data from chromosome 3 to 22. Narrowing down the parameter space that includes likely scenarios, we evaluated the goodness-of-fit of the models with the full genome data (i.e., chromosomes 1 to 22). Fitting of the model was measured as an approximate Bayes factor (aBF) (82) where the model with the highest likelihood was compared to each of all other models. We considered log₁₀-scaled aBF > 2.0 as decisive evidence on the data in favor of the highest likelihood model, as opposed to the other model (83).

To model admixture events in the prehistory of Japan, we applied qpWave v600 and qpAdm v1000 in the AdmixTools v6.0 package to the merged datasets of ancient East Eurasians and the SGDP populations (77, 84). Our analysis only used transversion sites with global minor allele frequencies of >1%, coupled with the option of “allsnps: YES.” We used a set of nine Eurasian populations, including Sardinian (n = 3), Kusunda (n = 2), Papuan (n = 14), Dai (n = 4), Ami (n = 2), Naxi (n = 3) (27), Tianyuan (n = 1) (85), Chokhopani (n = 1) (18), and Mal’ta (n = 1) (86), as an outgroup with or without a subset of Jomon individuals (JpKa6904, JpOd181, and IK002) in our modeling. We considered populations as sources of an admixture event only if the admixture between the sources is supported from modeling with and without a subset of Jomon in the right populations.

We used DATES v753 (45) to estimate the time of two different admixture events in the Kofun individuals: one is the admixture between the Jomon and Northeast Asian ancestry, and the other is between the Jomon and East Asian ancestry. For the former model, we used the West Liao River individuals labeled as WLR_BA_O and HMMH_MN as the second source. East Asian ancestry in the latter model was represented by Han Chinese in Beijing, China (CHB) in the 1000 Genome phase 3 data (28). The estimated date in generation was converted into years with the assumption of 25 years per generation, which was further added into a mean age of median dates from three Kofun individuals (i.e., 1348 years before present). The parameter settings that we used are as follows: binsize: 0.001, maxdis: 1.0, runmode: 1, mincount: 1, and lovalfit: 0.45. The SE was estimated from a weighted block jackknife method.