Ethics and inclusion statement

Permission to access the ancient DNA in the human remains in this study was approved by the archaeological team that lead the excavation from Shaanxi Academy of Archaeology, Institute of Archaeology (Chinese Academy of Social Sciences) and Archaeology Institute of National Museum of China. The Institutional Review Board of the Chinese Academy of Sciences, Institute of Vertebrate Paleontology and Paleoanthropology provided further monitoring and permission for the sampling of ancient humans in this research. All the work was done in collaboration with local archaeologists, who were named co-authors for their contributions to the collection of material and archaeological information, such as on-site photographs, classification of high-level tombs and/or discussions that contributed to the associations derived from the archaeological research cited in this study. All wet laboratory work and data analysis are performed with equipment from the Molecular Paleontology Laboratory, Institute of Vertebrate Paleontology and Paleoanthropology, Chinese Academy of Sciences.

Ancient DNA experiments and sequencing

We sampled and sequenced 207 human remains from Shaanxi and Shanxi provinces, China, among which 169 individuals were analysed in this study (Supplementary Table 1). All extraction, sequencing and data processing of ancient human samples were carried out in dedicated laboratories at the Institute of Vertebrate Paleontology and Paleoanthropology, Chinese Academy of Sciences, in Beijing. Following standard protocols⁴⁷, DNA was extracted from each sample from less than 100 mg of bone powder, obtained through drilling. We prepared double-stranded libraries (denoted ‘DS’ in Supplementary Table 1) for 134 samples with uracil-DNA glycosylase partially treated library protocol (denoted ‘half UDG’)^48,49 (Supplementary Table 1). For 35 samples, we prepared single-stranded libraries (denoted ‘SS’, Supplementary Table 1) with full UDG treatment (4 samples; denoted ‘UDG’, Supplementary Table 1) or no UDG treatment (31 samples; denoted ‘No’ in Supplementary Table 1). To collect enough DNA for capture, libraries were amplified for 35 cycles using the AccuPrime Pfx polymerase. We then evaluated the amount of DNA extracted per sample using a Thermo Scientific NanoDrop 2000 spectrometer. We applied a capture strategy on both mitochondrial and nuclear DNA. For mitochondrial DNA (mtDNA), we used oligonucleotide probes synthesized from the complete human mitochondrial genome⁵⁰, for nuclear DNA, oligonucleotide probes targeted 1.2 million SNPs (the ‘1,240k’ SNP panel) were applied. After enrichment, sequencing was performed on an Illumina MiSeq sequencing platform to generate 2 × 76 base pairs (bp) paired-end reads for the mtDNA and an Illumina Hiseq 4000 sequencing platform to generate 2 × 100 bp and 2 × 150 bp paired-end reads.

Read alignment and variant calling

We used leeHom⁵¹ to trim adaptors and merge paired-end reads into a single sequence (minimum overlap of 11 bp), keeping only merged reads with a length of at least 30 bp. Reads were aligned with BWA (v.0.5.10)⁵² using the bam2bam command with default parameters, except for samples with no UDG treatment, for which we used the parameters -n 0.01, -l 16500 and -o 2. We aligned the mtDNA reads to the revised Cambridge Reference Sequence⁵³ and the nuclear DNA reads to the human reference genome hg19 (ref. ⁵⁴). Duplicate reads with the same orientation, start and end positions were removed, and reads with a minimum mapping quality score of 30 were kept for analysis. The frequency of terminal C-to-T misincorporations was used to validate ancient DNA sequences, and contamination rates were estimated on the basis of two approaches. For all the individuals, we applied ContamMix⁵⁵ to compare mtDNA fragments between the new consensus mitochondrial genomes with the present-day sequences⁵⁰. To minimize the impact of damaged bases, we ignored the first and last five positions of the fragments during estimation. We treated the libraries as contaminated if the estimated contamination rate was greater than 5%. Contamination rates for men were also estimated using ANGSD⁵⁶, leveraging the fact that men have one copy of the X chromosome, and verified using HapCon⁵⁷, to improve the performance of low-coverage data. To keep enough individuals for further analysis, for 12 individuals with contamination above 5% (‘b’ annotated in the column of SNP number, Supplementary Table 1), we restricted our analysis to only damaged fragments with ancient DNA characteristics. The damaged fragments were obtained by pmdtools v.0.60 (ref. ⁵⁸) with the –customterminus parameter, keeping fragments with at least one C → T substitution in the first three positions at each end. To eliminate the potential bias caused by the terminal deaminated cytosines, we masked 2 bp at the end of mapped reads for all the double-strand libraries with half UDG treatment, and 5 bp at the end of the reads were masked for all the single-strand libraries with no UDG treatment. To generate pseudo-haploid genotypes, heterozygote SNPs were randomly sampled to determine a single allele for the individual. During genotyping, the first and last 5 or 2 positions of the fragments were ignored for non-UDG-treated and UDG-treated libraries, respectively, and 13 poorly covered samples (with fewer than 27,000 SNPs) were removed.

Uniparental haplogroup identification

Mitochondrial sequences for each individual were mapped to the revised Cambridge Reference Sequence⁵⁹. We only kept reads of a minimum of 30 bp in length and with a minimum mapping quality of 30. Haplogroups for each individual were called using HaploGrep2 (ref. ⁶⁰) based on PhyloTree Build v.17 (ref. ⁵⁹). We also confirmed all the haplogroups using the phylogenetic tree constructed with mtphyl v.5.003 and found that two individuals with an R# haplogroup (R + 16189)²⁷ were assigned into the subclade of Haplogroup B4c1a since the 9-bp deletion (8281–8289). In comparison, four individuals with B haplogroups were assigned to the ancestral haplogroup with the best fit (that is, B4a, B4b1 and B4c1b). For the male individuals, Y-chromosome haplogroups were determined by identifying the assigned position in the phylogenetic tree on the basis of the International Society of Genetic Genealogy dataset version 9.77 (www.isogg.org/tree). In cases in which the most derived allele upstream of the Y-chromosome was a C to T or G to A substitution, indicative of possible deamination, at least two derived alleles were required to assign the Y-chromosome haplogroup. Otherwise, the haplogroup of the tested individual would be assigned to the ancestral haplogroup. When the subclade of the haplogroup assignment could not be determined, the haplogroup of the individual would be assigned to the most recent ancestral haplogroup they best fit (for example, No).

Population structure analysis

We conducted a principal component analysis (PCA) with smartpca in the EIGENSOFT package⁶¹. To calculate the principal components, we used 82 present-day populations from the Affymetrix Human Origins dataset⁶². We merged newly sequenced and published ancient individuals to the Human Origins dataset and projected them using the following program settings: ‘lsqproject: YES’, numoutlieriter: 0, and shrinkmode: YES. Newly sequenced or previously published ancient individuals were projected onto the principal components calculated based on present-day Eurasians (Fig. 1c) or only the East Asians (Fig. 1d). We estimated individual ancestries by model-based maximum likelihood clustering using ADMIXTURE⁶³. We used 44 of the 82 populations used in the PCA, along with 10 present-day Han and Tibetan populations from ref. ⁶⁴. Before the admixture analysis, we pruned genotypes with high linkage disequilibrium (r² > 0.4) using PLINK (version v.1.90)⁶⁵ and the parameters ‘-indep-pairwise 200 25 0.4’ were applied for SNP filtering, leaving 597,573 SNPs. ADMIXTURE analysis was conducted with K from 2 to 10. For each K, we ran the analyses ten times with different seeds to estimate the cross-validation error, and the best K was determined according to the lowest cross-validation error.

Relatedness analyses

Kinship patterns among the samples from Shaanxi and Shanxi provinces were analysed using READ v.2 (ref. ⁶⁶) to determine the degree of kinship, and confirmed by lcMLkin⁶⁷. Further connections of the pedigrees or individuals were investigated using ancIBD⁶⁸. For Huangchengtai samples, we introduced the third hidden Markov-model-based approach by KIN⁶⁹ to test all the kinships within the second-degree. After filtering 22 genetically identical individuals, we found 25 kinship pairs with high confidence (Supplementary Table 3), consisting of 8 pairs of first-degree relationships, including 1 full sibling and 7 parent–offspring, and 17 second-degree relatives (Supplementary Table 3). For the READ analysis, a genome-wide approach was applied for calculating a single value (P0) across all sites without splitting the genome into windows, and the average P0 was then normalized by the median of all average pairwise P0 across all samples. To estimate standard errors, a block-jackknife approach was applied with blocks of 50 mega bp (Mbp). For the differentiation between parent–offspring and siblings, a different window size of 20 Mbp was used. We ran separate analyses by three groups: all Middle Neolithic Shaanxi samples, all Late Neolithic Shaanxi samples and all Shanxi samples. We used unrelated individuals without having first-degree and second-degree kinships estimated by READ for subsequent genetic analyses.

A further genotype likelihoods-based method to determine kinship, lcMLkin, was applied that considers the inaccuracy of genotype calling when sequence coverage is low. lcMLkin outputs the estimated probability of two diploid individuals sharing zero (k₀), one (k₁) or two (k₂) alleles that are identical by descent (IBD) and calculates the combined kinship coefficient by the equation: r = k₁/2 + k₂. The kinship categories (for example, identical twins or self, parent–offspring, full siblings, second-degree and unrelated) were determined by comparing with the theoretical expectation for k₀, k₁ and k₂ (ref. ⁷⁰). The method requires a SNP set with minor allele frequency higher than 5% and without linkage disequilibrium with each other. SNPs with allele frequency lower than 5% among all present-day East Asians from the Simons Genome Diversity Panel dataset⁷¹ were removed, and the resulting data were pruned for linkage disequilibrium using PLINK, with the parameter ‘-indep-pairwise 200 25 0.5’, which resulted in 135,642 SNPs available for the downstream analysis. We called genotype likelihoods at these SNP sites for our ancient individuals using the script SNPbam2vcf.py available with lcMLkin and estimated their biological relatedness using lcMLkin (Supplementary Table 3). When two samples from the same burial site were identified as identical twins or the same individual, the sample with the lower coverage was removed from analysis. For samples with first-degree or second-degree familial relationship, we only retained genome-wide data for the individual who had the higher coverage and no first-degree or second-degree relationships with the remaining individuals in the same site for further downstream genetic analysis. This resulted in 25 samples being excluded from the population analyses.

For the shared IBD analysis, genomes from individuals with coverage above 0.5× were imputed by GLIMPSE2 (refs. ^72,73) setting parameters for quality control in phasing as –mapq 30 and –baseq 30. The 1000 Genomes Phase 3 dataset was used as a reference panel and was downloaded at http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/. To obtain the IBD sharing information, imputed files were merged by chromosomes and analysed by the software ancIBD⁶⁸ with suggested parameters, including genotype posterior probabilities higher than 0.99, gap maximum at 0.0075 and IBD blocks with more than 220 SNPs. The results of IBD sharing between pairs of individuals are recorded in Supplementary Table 12. We focused on IBD connections between individual pairs both having coverages above 1×. We also provided the results for IBD connections when the coverage of either of the individuals in a pair is between 0.5 and 1×, but denoted them as ‘low confidence’. To confirm the within-second-degree kinship relationships estimated by READ, lcMLkin and KIN (Supplementary Table 3 and Supplementary Figs. 16 and 17a), we counted and plotted the sum and number of IBD segments longer than 12 cM, which could be applied to infer kinship relationships⁶⁸. We also depicted IBD segments that were longer than 8 cM (Supplementary Figs. 17b–d) in karyotype plots for these individuals who shared within-second-degree kinships.

To further infer distant kinships of more than second degrees (denoted as third-degree to fifth-degree degree kinships in Fig. 3, Extended Data Fig. 3 and Supplementary Figs. 18–21), we counted and plotted the length distribution and the karyotype of shared IBDs longer than 8 cM. The kinship relations were determined based on the goodness of fit between the observed and the expected curves of various kinship categories, denoted in the top right of each karyotype plot (Supplementary Figs. 17–21). Individual pairs within each site with IBD sharing shown in figures were denoted as possibly having third-degree to fifth-degree kinship relationships (Figs. 3 and 4 and Supplementary Figs. 17–21). Individual pairs across sites sharing single IBDs of at least 25 cM in length were also counted (Extended Data Fig. 3). We also compared the sum and number of IBDs that the segments longer than 8 cM or 12 cM between individuals with different social status (that is, tomb owners and sacrificed victims) at different sites (Supplementary Table 13 and Extended Data Fig. 3). We only included individuals with an average SNP panel coverage of at least 1×. Only one individual was considered for pairs identified as genetically identical.

ROH

The presence of close-kin mating and kinship-based mating systems could be reflected by the ROH, which represent long stretches of homozygous segments along the genome of an individual. For 171 individuals with SNP counts at least 29,000, we applied hapROH⁷⁴ that detects ROH in low-coverage ancient DNA data using haplotype information from a modern phased reference panel. We detected ROH with four length classes (4–8 cM, 8–12 cM, 12–20 cM and more than 20 cM). Results for individuals with SNP counts more than 400,000 were recorded as ‘high confidence’, whereas others were recorded as low confidence (Supplementary Table 14). Observed and expected ROH blocks above 4 cM were plotted using hapROH.

Genetic clustering among new samples

Samples were grouped using a combination of D statistics and PCA. We calculated the D statistics in form D(Sample 1, Sample 2; Population, Mbuti) for each pair of samples. Here, we took Mbuti as the outgroup. Sample 1 or sample 2 are two unrelated individuals from each archaeological site. The 45 individuals and populations (both ancient and present-day) making up the population in the above statistic are grouped as follows:

aNorthEA(14)

AR_EN, Bianbian, DevilsCave_N, Chokhopani, HMMH_MN, Kolyma, Mebrak, Mongolia_N_East, Okunevo_EMBA, Shamanka_EN, WLR_MN, WLR_LN, WLR_BA, Yumin; aSouthEA(7): Liangdao1, Liangdao2, Longlin, Man_Bac, Nui_Nap, Qihe and Qihe3.

DeepAsia(9)

G1, Jomon, Malta1, Onge, Papuan, Tianyuan, USR1, Vanuatu, Yana; aWest/SouthAsia(15): Afanasievo, Anatolia_N, Botai_CA, Ganj_Dareh_N, IndusPeriphery, Hajji_Firuz_C, Harappan, Iran_N, Karelia, Kotias, Kostenki14, Shahr_I_Sokhta_BA3, Ust-Ishim, Vestonice16 and Yamnaya_Kalmykia.SG.

If the samples can be grouped in the same genetic cluster, we predict that D ≈ 0 (|Z score| < 3) for most populations. Samples that deviated from this expectation within groups, or were found to be outliers in the PCA, were separated from the main group. Results were only considered from sample pairs having at least 25,000 overlapping SNPs. We summarize the matrix for the count of significant pairwise D statistics for each archaeological site in Supplementary Table 2. The detailed genetic grouping is listed in Supplementary Table 1 and the PCA clustering is presented in Fig. 1. The 18 populations defined in the text are denoted and highlighted in grey in Supplementary Table 4.

Outgroup f
₃ and D statistics

We calculated f₃ statistics using qp3Pop (v.412) with the form f₃(Population X, Population Y; Mbuti), measuring the shared genetic drift between all combinations of populations relative to the outgroup. We used the present-day central African population, Mbuti, as the outgroup and compared newly sequenced and previously published ancient populations within or outside East Asia (Extended Data Fig. 1). The higher the f₃ statistic, the more genetic drift (or shared genetic similarity) two populations share relative to Mbuti. We calculated D statistics using qpDstat (version 712) with the form D(population X, population Y; population Z, Mbuti), measuring the shared number of alleles between all combinations of grouped new populations and a diverse array of previously published ancient and present-day populations (Supplementary Tables 4, 5 and 7–11). A negative D statistic means that population Y shares more alleles with population Z than it does with population X. A positive D statistic means that population X shares more alleles with population Z than it does with population Y. Both the outgroup f₃ and D statistics were calculated using AdmixTools⁷⁵.

Admixture modelling with qpAdm

The ancestry proportions of ancient populations were estimated using qpAdm (v.634) in the AdmixTools package, modelling one, two or three different sources. Distal and proximal modelling were used to model the ancestry of target populations, in which two modelling types differ in the relative age of the source populations. Distal modelling considers older source populations with larger genetic distance (Yumin, sEastAsia_EN, Xitoucun, Man_Bac, Coastal_nEastAsia_EN, Xingyi_EN, AR14K, YR_MN, Wuzhuangguoliang and WLR_MN) and proximal modelling looks at younger source populations with closer genetic distances (Yumin, Wuzhuangguoliang and YR_MN), which was applied for the outlier modelling (sEastAsian_EN = Qihe2 and Liangdao2; Coastal_nEastAsia_EN = Bianbian, Boshan and Xiaogao)³². Both model types used the same set of reference populations, which include Mota, Ust-Ishim, Kostenki14, Iran_N, IndusPeriphery, LBK_EN, Motala12, Kotias, AR33K, Yana, Karelia (IndusPeriphery = Shahr_I_Sokhta_BA2 and Shahr_I_Sokhta_BA3 (5–4 ka) from Shahr-i-Sokhta in Iran, and Gonur Depe (Gonur2_BA) (roughly 4 ka) from the Bactria-Margiana Archaeological Complex in Turkmenistan)^32,76. As for the subgroups (Shimao_HJGD3 and Zhoujiazhuang3) potentially having connections with Western or Central Asian populations, we consider extra sources (Anatolia_N, Afanasievo, AYTH, Saidu_Sharif_H, Kashkarchi_BA, Zevakinskiy_LBA, IndusPeriphery, Bustan_BA, Botai_CA, Satsurblia, Gonur1_BA, Dzharkutan1_BA, Gogdara_IA, Loebanr_IA, Shahr_I_Sokhta_BA1, LaBrana1, Levant_N, Stuttgart and Loschbour) along with the source populations in the distal model described above, and with the same reference population but excluding Iran_N and IndusPeriphery. A ‘rotating’ scheme⁷⁷ of source and reference populations was used for distal and proximal modelling. Wuzhuangguoliang and the extra sources were used only as the source population. As the genetic make-up of the proximal sources is too close to compete for the most optimal model, on the basis of the determined sources (Wuzhuangguoliang or further sources) from distal modelling, we fixed them as the source populations. The other parameters we adopted for qpAdm modelling are: ‘allsnps: YES’, ‘details: YES’ and ‘summary: YES’. The model was deemed plausible if the tail probability of rank0 is above 0.05 and the estimated admixture proportions are between 0 and 1. A valid model with more source populations was considered only when fewer sources were rejected, and results with the lower number of sources were marked as the high-confidence model in Supplementary Table 6.

We observed two-way admixture models for Shimao_HJGD3, including 87–93% Yellow River Yangshao ancestry (represented by Wuzhuangguoliang or YR_MN) with an extra 7–13% ancestry components from diverse non-East Asian ancestries, represented mainly by Iranian farmer-related and Western Hunter-gatherer-related ancestry. Zhoujiahzuang3 could still be modelled as having one source from either Coastal_nEastAsia_EN or Wuzhuangguoliang (Supplementary Table 6). To get a broader source for the Wuzhuangguoliang population, we also conducted further qpAdm testing by swapping AR14K with the later ARpost9K population, and we found ARpost9K could better represent the AR14K and Yumin ancestry components for Wuzhuangguoliang ancestry.

Modelling of Wuzhuangguoliang ancestry

To further understand the composition of ancestry sources for Wuzhuangguoliang, we simulated the admixed populations with 20 replicates to better represent the potential ancestries with admixture proportions from 0% to 100% in increments of 1% following the method implemented in ref. ⁷⁸. The measurement is based on f₄ statistics in the form of f₄(Wuzhuangguoliang, simulated mixture populations; nEA/outEA/steppe/AR-related, Mbuti) with a two-sourced mixture of Early Neolithic Shandong including Bianbian, Boshan and Xiaogao subpopulations with YR farming ancestry (YR_MN), or a three-sourced mixture by adding the third group with Yumin or Amur River ancestry (AR14K and ARpost9K) in comparison to nEA (Miaozigou_MN, WLR_MN and AR-related such as AR14K, AR19K, ARpost9K and DevilsCave_N), outside East Asian (outEA = Shamanka_EN, Lokomotiv_EN and Loschbour), steppe-related (Mongolia_N_East, Mongolia_N_North and Yumin) and Tibetans (Shannan2k and Zongri5.1k). The ancestral component selection is mainly based on the valid models from the qpAdm analysis. For the proportional simulation of the ternary simulated population, we first use the admixed population from YR and Coastal_nEA_EN/ARpost9K as the first admixed component, varying proportionally from 0% to 100%, and then add the second component of ARpost9K, Yumin or WLR_MN proportionally from 0% to 100%. The ranges of proportion thresholds that marked a plausible range of the second or third admixed ancestry are estimated by qpAdm modelling (Supplementary Table 6) with the source proportion ± standard error rate. The past models of admixed ancestries of Wuzhuangguoliang could be explained by Wuzhuangguoliang harbouring an ancestry that contains more affinity with Neolithic Shandong, Amur River, West Liao River or Tibetan ancestries. Although none of the specific nEA groups served as a good proxy for this unknown ancestry because adding them as the second or third source was still insufficient to model Wuzhuangguoliang ancestry (Supplementary Figs. 12–15).

Demographic modelling with qpGraph

We applied qpGraph and findGraphs functions from AdmixTools and AdmixTools2 packages⁷⁹, respectively. First, we followed a general approach for modelling admixture graphs using qpGraph (v.6065) in the AdmixTools package, which started with a basic and well-understood tree (including the central African Mbuti as an outgroup, the early western Eurasian UstIshi and early East Asian Tianyuan). We then added extra populations (sEastAsia_EN, YR_MN, Wuzhuangguoliang, Shimao_HJGD1, Yumin, Wuzhuangguoliang_o and Xinhua_o) one at a time in their best-fitting positions iteratively³². An optimum tree model could be constructed based on the observed f statistics (f₂, f₃ and f₄ for all possible pairs of populations). We required a |Z| score of less than 3 between the observed and expected values (determined by the block Jackknife) to accept the model. A small number (0.0001) was added to the diagonal entries of the estimated covariance matrix of the f statistics (Q matrix) to stabilize the matrix inversion. The qpGraph program was run with the following recommended parameters⁸⁰: ‘outpop: Mbuti, blgsize: 0.05, lsqmode: YES, diag: 0.0001, hires: YES, initmix: 1000, precision: 0.0001, zthresh: 0, terse: NO, useallsnps: NO’. To not miss any unexplored graph models, we first carried out fully automated graph exploration using findGraphs tools, allowing 0 to 8 admixture events to occur in 100 algorithm iterations (Stop_gen = 100). Next, we constrained deeply diverged populations (Mbuti, Tianyuan, sEastAsia_EN, Ust-Ishim), assuming they are non-admixed with 100 algorithm iterations, setting admixture events from 0 to 3. Here sEastAsia_EN, including Qihe2 and Liangdao2, and all the denoted populations are similar to the admixture graph by qpGraph in AdmixTools. Graphs of best-fitting models are listed in Fig. 2 and Supplementary Figs. 8 and 9.

Treemix analysis

The phylogenetic relationships among populations were estimated using Treemix v.1.13 (ref. ⁸¹). The dataset for the Treemix analysis included the newly sequenced cohort including five groups based on their genetic and spatiotemporal characteristics: preShimao_5k, Shimao_4k, Taosi_4k, Shimao_4k_o and preShimao_5k_o (denoted in Supplementary Table 1) and previously published cohorts including DevilsCave_N, WLR_MN, WLR_LN, Coastal_nEastAsia_EN, Yumin, AR19K, sEastAsia_LN, sEastAsia_EN, Tianyuan, Yana and Ustlshim. The maximum likelihood tree was rooted by Mbuti (-root Mbuti) and linkage disequilibrium was compensated for by grouping sites in blocks of 500 SNPs (-k 500). A round of global rearrangements and no sample size correction were used with the parameters ‘-global -noss’, allowing 0 to 7 migration events (-m 0–7). For m = 2, we ran 1,000 bootstraps for the maximum likelihood tree (-bootstrap), then assessed 1,000 bootstrapped trees in phylip, using the consense command⁸² to count the number of times a particular branch or populations clustered for maximum likelihood tree with migration. Results are shown in Fig. 2. The inferred maximum likelihood trees with migrations from 0 to 7 and the corresponding residuals (Supplementary Fig. 6) were visualized with an R script from Treemix v.1.13.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.