Jomon: Paleolithic Contingent in Modern Japanese

new

Admin
Administrator

Posts: 81,722

Jomon: Paleolithic Contingent in Modern Japanese Sept 21, 2021 22:40:41 GMT

Quote

Post by Admin on Sept 21, 2021 22:40:41 GMT

Materials and Methods
Samples
The skeletons that were analyzed in the current study were excavated in 1998 (Nishimoto, 2000) from the Funadomari shell mound, which is located in the northern part of Rebun Island (Figure 1). The two skeletons were housed at Sapporo Medical University and were dated to 3500–3800 years ago based on C14 dating of charcoal from the hearth (Nishimoto, 2000). From the stable isotope value for one of the samples (F23), we inferred the source of diet to be a mixture of marine and land animals. We calculated F marine from the proportion of marine animal (δ 13C = −12.5‰) and land animal (δ 13C = −21.0‰). The age of F23 was also dated to be 3550–3960 cal BP (Paleo Labo Co., Ltd.) (Supplementary Table 1). These values are consistent with “the late Jomon period” according to pottery-based dating (Nishimoto, 2000). We extracted DNA from the left second and right first mandible molar teeth of Funadomari No. 5 (F5) and No. 23 (F23), respectively. Their DNA was exceptionally well preserved (Adachi et al., 2009). In the previous study (Adachi et al., 2009), we used the dental pulp and the dentin around the dental cavity for DNA extraction. These DNA samples were obtained by drilling the tooth of each skeleton from the apical foramen. In the present study, we used the residual tooth powder of F5 and F23 that have been stored at −70°C at the Department of Legal Medicine, University of Yamanashi since 2009.

Figure 1
Map of East Asia. (a) Geographical location of Rebun Island, Hokkaido, Japan. (b) Geographical location of Funadomari site.

DNA extraction and library preparation
DNA from F5 and F23 (Figure 2) was extracted according to the methods described by Adachi et al. (2018) in the laboratory dedicated specifically for ancient DNA analyses at the University of Yamanashi. Briefly, the powdered sample was decalcified with 0.5 M EDTA (pH 8.0) at 56°C for 24 h, and the buffer was replaced with a fresh solution and decalcified again at 56°C for 36 h. Thereafter, the decalcified samples were lysed in 1000 μl of Genomic Lyse buffer (Genetic ID, Fairfield, IA, USA) and digested with 50 μl of 20 mg/ml proteinase K at 56°C overnight. Finally, the lysate was treated twice with 1500 μl phenol–chloroform–isoamyl alcohol (25:24:1) (Thermo Fisher Scientific, Waltham, MA, USA) and 1500 μl chloroform (Wako Pure Chemical Industries, Osaka, Japan). DNA was extracted from the lysate using a FAST ID DNA extraction kit (Genetic ID) in accordance with the manufacturer’s protocol.

Figure 2
Photo of female 23 (F23), whose genomic DNA sequence was determined. Scale bar = 10 cm.

Following the protocol of Shinoda et al. (2017), which was modified from the protocol of Kanzawa-Kiriyama et al. (2017), three and seven DNA libraries were prepared from DNA extracts of F5 and F23, respectively. For creating blunt-ends and adding a terminal 3′ adenine and adapter ligations, we used a GS Titanium Rapid Library Preparation kit (454 Life Science Corporation, Branford, CT, USA). The conditions for the generation of blunt-ends and addition of adenine were 20 min at 25°C, 20 min at 12°C, and 20 min at 72°C, with a hold at 4°C. For the adapter ligation, we used 1 μl of Illumina LT Set A adapter mix diluted to 10× using water or 2 μl of 10 μM short adapter designed by Shinoda et al. (2017) with 7 μl of 50% PEG-4000 (Hampton Research, Aliso Viejo, CA, USA). The short adapter was prepared by hybridizing two oligonucleotides (Short_Adapter_1: 5′-ACACTCTTTCCCTACACGACGCTCTTCCGATC-phosphorothioate-T-3′ and Short_Adapter_2: 5′-PHO-GATCGGAAGAGCACACGTCTGAACTCCAGTCAC-3′) as follows. In a polymerase chain reaction (PCR) tube, 20 μl of 500 μM CL53, 20 μl of 500 μM CL73, 9.5 μl TE (10 mM Tris–HCl, pH 8.0), and 0.5 μl of 5 M NaCl were combined. This mixture was incubated for 10 s at 95°C in a thermal cycler and cooled to 14°C at 0.1°C/s. A final concentration of 100 μM was obtained by dilution with 50 μl TE. The ligation was carried out for 3 h at 25°C. We purified the post-ligated solutions using the QIAquick PCR purification kit (Qiagen, Hilden, Germany) and eluted the ligated samples in 30 μl elution buffer (EB). When using the short-adapter, we diluted the post-ligated solution 3-fold with water and then used AMpureXP (Beckman Coulter, Brea, CA, USA) for purification. The entire volume of the purified library was used for PCR to complete the library preparation. We mixed 50 μl of the DNA library, 50 μl Multiplex PCR Master Mix (Qiagen), and 1000 nM of each multiplexing PCR primer and dispensed the mixture into two 200 μl single PCR tubes. Cycling conditions were as follows: 15 min at 95°C; 12 cycles of 30 s at 95°C, 30 s at 65°C, and 30 s at 72°C, with a final extension at 72°C for 10 min. The first PCR products were purified with an AMpureXP kit and eluted in 50 μl EB. A second-round of PCR was set up as follows: 10 μl product from the first PCR, 10 μl of 10× AccuPrime Pfx Reaction Mix, 1000 nM of each multiplexing PCR primer, 2.5 U of AccuPrime Pfx DNA polymerase (Thermo Fisher Scientific), and H2O were mixed to obtain a final reaction volume of 100 μl. We dispensed the mixture into two 200 μl PCR tubes. Cycling conditions were as follows: 95°C for 2 min; 8–12 cycles of 95°C for 15 s, 60°C for 30 s, and 68°C for 30 s; and a hold at 4°C. The second PCR products were purified with a QIAquick PCR purification kit and eluted in 50 μl EB. The purified libraries were quantified with an Agilent 2100 Bioanalyzer DNA High-Sensitivity chip (Agilent Technologies, Santa Clara, CA, USA) for sequencing.

Sequencing and sequence mapping
DNA sequencing was conducted for the DNA libraries on an Illumina 2500 platform using a paired-end run with 100 cycles at Hokkaido System Science Co., Ltd. (Sapporo, Japan). After trimming the adapter sequence with CutAdapt (Martin, 2011) (−q 25 −m 35 -a AGATCGGAAGAGC -A AGATCGGAAGAGC), the remaining paired reads were mapped to the human reference genome (hs37d5) using BWA (http://bio-bwa.sourceforge.net/index.shtml; mem -M) (Li and Durbin, 2009). The frequencies of mapped reads were examined using SAMtools (http://samtools.sourceforge.net/index.shtml), with the flagstat option (Li et al., 2009). Mapped reads with a mapping quality of 30 or more were retained, sorted, and merged. We then checked the insert size of paired reads, and exceptionally long reads were removed to exclude mismapped reads (≤400 for F5 and ≤500 for F23). PCR duplicates were removed using Picard MarkDuplicates (version 1.119; broadinstitute.github.io/picard/). After eliminating hard-clipped and soft-clipped reads and trimming three bases in the termini of the reads using BamUtil software (http://genome.sph.umich.edu/wiki/BamUtil), the resulting aligned reads were realigned using GATK IndelRealigner (--consensusDeterminationModel KNOWNS_ONLY -LOD 0.4) and recalibrated using GATK BaseRecalibrator. We then examined read length distributions for F5 and F23 based on the insert size.

Metagenomic analysis
We performed metagenomic analysis to check for nonhuman DNA contamination. The reads were first preprocessed with the ‘tantan’ repeat masker (Frith, 2011) using options ‘x N –r 0.0005.’ The reads were then mapped to the reference human genome build GRCh38.p11, accession GCF_000001405.37, from RefSeq (O’Leary et al., 2016) with minimap2 (Li, 2018) using options ‘–x sr.’ The reads that failed to map were extracted as the ‘nonhuman’ read set. A subset of 100000 of such reads (about 10 Mbp) from each sample was then analyzed with BLAST (https://www.ncbi.nlm.nih.gov/BLAST/) using parameters ‘blastn –task blastn–dust no –evalue 3.80e–2 –dbsize 3200000000’ (Altschul et al., 1997). The BLAST search was performed against available GenBank (NCBI Resource Coordinators, 2018) and RefSeq (O’Leary et al., 2016) genome assemblies, obtained via the GenomeSync database (http://genomesync.org/). Overall, 143139 genomes were used, totaling about 1.73 Tbp of sequence. The search was performed and analyzed with the Genome Search Toolkit (http://kirill-kryukov.com/study/tools/gstk/). The reads were taxonomically classified by determining their specificity taxa; a read was labeled as specific to a taxon if it had more significant BLAST hits to genomes within this taxon compared with its hits in all other genomes. KronaTools (Ondov et al., 2011) was used for visualizing the results.

Mitochondrial genome analysis
We used GATK HaplotypeCaller to call variant sites for mitochondrial genomes (-stand_call_conf 30 -stand_emit_conf 10). We manually compared the list of variants with SNPs reported in phylotree-Build 17 (van Oven and Kayser, 2009) to determine mtDNA haplogroups. Based on the haplogroups, the resulting list of variants was double-checked using IGV software (Thorvaldsdóttir et al., 2013). Mitochondrial genome consensus sequences were prepared using GATK FastaAlternateReferenceMaker and used for later contamination estimates.

Sex determination and Y chromosome analysis
We first determined the sex of F5 and F23 by counting the number of sequence reads mapped to X and Y chromosomes. The size of the X chromosome is 10 times larger than that of the Y chromosome in the human reference genome, meaning that if the ancient individual was male, the X:Y ratio of mapped reads would be close to 10:1.

For Y chromosome analysis, we retrieved the human Y chromosome nucleotide sites registered at the International Society of Genetic Genealogy (ISOGG) website (https://isogg.org/tree/ISOGG_YDNA_SNP_Index.html) and then prepared an mpileup file of F5 overlapping with the listed nucleotide sites. We then determined the haplogroup classification of F5 and compared the F5 Y chromosome sequence and 60505 SNP sites listed in 1000 genome project Y chromosome data (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chrY.phase3_integrated_v2a.20130502.genotypes.vcf.gz); 24391 SNP site data were shared between them. A phylogenetic tree was generated for these sequences using MEGA7 (Kumar et al., 2016) with the following parameters: Tamura-Nei model, uniform rates, using all sites.

DNA authenticity analysis and error rate estimation
To clarify the authenticity of the mapped sequence reads, we checked the frequency of postmortem misincorporation and depurination, which are characteristic of ancient DNA, with mapDamage2.0 (Jónsson et al., 2013). After confirming the authenticity, we estimated the contamination frequency for mtDNA and for nuclear DNA sequences. For mtDNA, we followed the methods used by Kanzawa-Kiriyama et al. (2017), in which the consistency of reads mapping to the SNPs of mtDNA haplogroups and individual specific mutations was estimated. We also used Schmutzi software (cont-Deam.pl --library single) (Renaud et al., 2015), a tool that uses a Bayesian maximum as an a posteriori algorithm to estimate mtDNA contamination. For nuclear DNA, we used ANGSD software (‘MoM’ estimate from ‘Method 1’ and the ‘new_llh’ likelihood computation) (Korneliussen et al., 2014), which estimates contamination based on the X chromosome polymorphism rate in males. We followed the criteria of Fu et al. (2016), stating that libraries are effectively uncontaminated if their X chromosome contamination estimates are less than 2.5%. We also estimated the error rate based on the excess of derived alleles in ancient samples using ANGSD software (-doAncError -anc hg19ancNoChr.fa -ref hs37d5.fa). For error rate estimates, data in which 3–10 bases were trimmed from sequence termini were used.

Admin
Administrator

Posts: 81,722

Jomon: Paleolithic Contingent in Modern Japanese Sept 22, 2021 4:55:32 GMT

Quote

Post by Admin on Sept 22, 2021 4:55:32 GMT

Haplotype-based analysis
Panel construction
We estimated the haplotype phase of F23 with a statistical approach. We first downloaded seven Tibetan (of seven regional Tibetan populations, we randomly sampled one individual per population) and five Sherpa genomes from the Genome Sequence Archive (GSA: PRJCA000246) (Lu et al., 2016) and then mapped and filtered the sequences by following the workflow of Mallick et al. (2016). To obtain an accurate result, we phased the F23 genotype data with diverse population genotype data. Thus, genotypes were determined for F23 in addition to 278 samples from the Simons Genome Diversity Project (SGDP) (Mallick et al., 2016), seven and five samples of Tibetan and Sherpa, respectively (Lu et al., 2016), Altai Neanderthal (Prüfer et al., 2014), and Denisovan (Meyer et al., 2012) using the HaplotypeCaller algorithm implemented in GATK3.8 with the ‘-ERC GVCF’ option. The second method for genotype calling was conducted with gVCF format files using the GenotypeGVCF program of GATK3.8. We then applied the variant quality score recalibration (VQSR) approach to obtain the same degree of genotype call sensitivity between F23 and modern humans. VQSR was applied to the raw genotype call set using the Hapmap and 1000 Genomes Omni2.5 sites as true datasets, high-confidence SNPs from 1000 genome sites as training datasets, and dbSNP138 sites as known datasets. The variant sites were filtered at 99% sensitivity. In addition, for Altai Neanderthal and Denisovan, genotypes with genotype quality (GQ) less than 20 were masked as no call variants. We further removed the variant sites for haplotype phasing when the call rate was less than 97%, the genotype frequency deviated from the Hardy–Weinberg equilibrium (P < 0.00001), or only one heterozygous genotype was observed among the whole data set.

Identical by descent (IBD) and homozygosity by descent (HBD) detection
The segments of IBD and HBD were detected using IBD-seq (Browning and Browning, 2013) with default settings (ibdlod = 3.0, r2window = 500, and r2max = 0.15). The physical lengths of each IBD and HBD segment (in base pairs) were converted into genetic distances (in centimorgans, cM) by extrapolating the genetic map provided by HapMap (International HapMap Consortium, 2005). The lengths of IBD segments for each sample were summarized for each paired sample. These procedures were also applied to the simulated data set to see the impact of small population size on HBD length distribution. Thus, we simulated a scenario of a population with a small population size after Out-of-Africa. The population size was 20000 before 7700 generations ago and then decreased to 10000, 5000, and 200 at 7700, 2500, and 100 generations ago, respectively. The coalescent simulation was conducted with MaCS (Chen et al., 2009) to generate 10 samples (20 haploids). The genome was assumed to have the same length as hg19, and the recombination map of HapMap was used. The mutation rate was assumed to be constant across the genome (10−8/bp/generation). Differences in the distributions of IBD segments between two samples were tested by the two-sample Kolmogorov–Smirnov test using the ks.test function implemented in R version 3.4.2.

Human leukocyte antigen (HLA) allele inference
Read sequences were extracted from BAM files of F23 together with unmapped reads. Using the resulting FASTQ files, HLA alleles were estimated with HLA-VBseq (Nariai et al., 2015) with reference sequences provided by the IMGT/HLA database (release 3.29). We selected HLA alleles with the top two highest read coverages as the candidate HLA alleles for each locus. Alignment results of the HLA genes were reconfirmed by de novo assembly. The sequence reads mapped on the HLA genes were extracted from bam files in fastq format. The HLA-filtered reads were realigned to the target HLA genes of HLA-A, -B and -C, and allelic sequences of HLA genes were phased by de novo assembly using Sequencher DNA Sequence Analysis Software version 5.0.1 (Gene Codes Corporation, MI, USA). HLA genotyping was performed by BLAST search with IMGT/HLA database. The misalignments among similar sequence reads of HLA genes were realigned for specific HLA allele sequences. The aligned reads were manually examined with IGV software to double-check the candidate HLA alleles. Allele frequencies of modern populations were obtained from the Allele Frequency Net Database (http://www.allelefrequencies.net/default.asp).

Genotype calling of F23
We used GATK UnifiedGenotyper to call genotypes for all genomic sites (-dcov 600 -stand_call_conf 30 or 50-stand_emit_conf 10 --output_mode EMIT_ALL_SITES--genotype_likelihoods_model BOTH). Next, we masked NUMT (Shinoda et al., 2017) and tandem repeat regions (Prüfer et al., 2014) to restrict data analysis to unique regions. Then, we checked the distribution of depth of coverage to determine whether the depth was sufficient for genotyping (≥10). To investigate the effects of residual deamination in the genotyping quality, we checked the allelic distribution in heterozygous sites. If each allele was sequenced randomly at each heterozygous site, we expected a binomial distribution centered in the middle.

Lazaridis et al. (2014) pointed out that, for sites likely to be heterozygous, two additional peaks would be observed at low allele frequencies due to sequencing errors, contamination, and deamination in some ancient samples, which are originally homogeneous. Therefore, for heterozygote sites, -stand_call_conf 50 was used, and allelic balances greater than 0.2 were removed (Rasmussen et al., 2014). We also expected that segmental duplications and copy number variants (CNVs) would cause the depth and heterozygote frequency to be artificially high. We removed sites with a depth higher than 79 (2.13% of total covered sites) and lower than 10 (4.00% of total covered sites). In total, 6.13% of covered sites were removed.

Phenotyping of anthropological traits and detection of pathogenic mutations
Phenotypes of F23 for some anthropological traits (alcohol tolerance, earwax type, hair thickness, shovel-shaped incisors, iris color, and skin reflectance) were inferred from genotypes of the associated SNPs, as reported in previous studies (Harada et al., 1981; Edenburg and Bosron, 1997; Yoshiura et al., 2006; Motokawa et al., 2007, 2008; Fujimoto, 2008; Kimura et al., 2009; Yamaguchi et al., 2012). The height tendency of the F23 female was also estimated from genotypes of SNPs reported as genetic factors for adult height in a recent meta-analysis in East Asians (He et al., 2015). Seventy-five SNPs with P < 1 × 10−5 in stage 1 and combined P < 5.0 × 10−8 were used for height estimation. We used effect size in stage 2 because of the smaller SE and assumed that a person without any positive and negative effect alleles was of average height (Z = 0). In addition, ABO and Rh blood types were estimated using BOOGIE (Giollo et al., 2015).

Clinically significant variants were detected using Clin-Var entries based on the ANNOVAR annotation tool (Landrum et al., 2016; Wang et al., 2010). For detection of pathogenic rare amino acid substitutions, we used ClinVar to obtain information about relationships among variations and human health and data for 2049 Japanese individuals from the Integrative Japanese Genome Variation Database (Yamaguchi-Kabata et al., 2015) (https://ijgvd.megabank.tohoku.ac.jp/about/) and for 1208 Japanese individuals from the Human Genetic Variation Database (Higasa et al., 2016) (http://www.genome.med.kyoto-u.ac.jp/SnpDB/index.html), Exome Aggregation Consortium version 0.3.1 (Lek et al., 2016) (http://exac.broadinstitute.org), and 1000 Genomes Project database Phase 3 (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/) to narrow down rare variants. The functional effects of the variants were evaluated using combined annotation-dependent depletion (CADD) for scoring the deleteriousness of non-coding variants as the clinical validity (Kircher et al., 2014) and Protein Variation Effect Analyzer (PROVEAN) to predict whether an amino acid substitution or indel had an impact on the biological function of a protein (Choi et al., 2012).

Data preparation for modern and other ancient human genomes
Seven data sets from the 1000 Genomes Project Consortium (2012), the SGDP (Mallick et al., 2016), genomes of Tibetan and Sherpa (Lu et al., 2016), genotypes of fully public Affymetrix Human Origins present-day individuals and ancient individuals analyzed by Lazaridis et al. (2016) and Skoglund et al. (2016) (the Pacific dataset), genotypes of three modern Japanese populations (Japanese Archipelago Human Population Genetics Consortium, 2012), and data from the HUGO Pan-Asian SNP Consortium et al. (2009) were used (Table 1). For filtering the SGDP dataset, we used SGDP-lite version 3 (http://sharehost.hms.harvard.edu/genetics/reich_lab/cteam_lite_public3.tar), for which nine filter levels were assigned for the variant call. We chose filter level 1, which is believed to be useful for most applications and achieves a good balance between sensitivity and low error rate and retains only the transversion sites. We used cTools (https://github.com/DReichLab/cTools>) for the filtering, retained only variant sites and merged the data with 1000 Genomes Project data. We used LiftOver (http://hgdownload.cse.ucsc.edu/admin/exe/macOSX.x86_64/liftOver) to convert hg18 coordinates to hg19 coordinates for the data from the Japanese Archipelago Human Population Genetics Consortium (2012).

Admin
Administrator

Posts: 81,722

Jomon: Paleolithic Contingent in Modern Japanese Sept 22, 2021 22:00:17 GMT

Quote

Post by Admin on Sept 22, 2021 22:00:17 GMT


Table 1 List of datasets and populations used in the current study
Dataset	Array or whole genome	Regions or contries	Populations	Reference
Funadomari Jomon	Genome	Japan	3500–3800 YBP Hokkaido Jomon (DRA008031)	This paper
Sanganji Jomon	Genome	Japan	3000 YBP Tohoku Jomon	Kanzawa-Kiriyama et al. (2017)
Mal’ta MA1	Genome	South central Siberia	24000 YBP south-central Siberian	Raghavan et al. (2013)
Tianyuan	Genome	East Asia	40000 YBP East Asian	Yang et al. (2017)
1000 Genomes	Genome	Africa	LWK (Luhya in Webuye, Kenya)	The 1000 Genomes Project Consortium (2012)
Europe	CEU (Utah Residents (CEPH) with Northern and Western European Ancestry)	
East Asia	JPT (Japanese in Tokyo, Japan), CHB (Han Chinese in Beijing, China), CHS (Southern Han Chinese), CDX (Chinese Dai in Xishuangbanna, China), KHV (Kinh in Ho Chi Minh City, Vietnam)	
SGDP (Simons Genome Diversity Project)	Genome	Africa	BantuHerero, BantuKenya, BantuTswana, Biaka, Dinka, Esan, Gambian, Ju hoan North, Khomani San, Luhya, Luo, Mandenka, Masai, Mbuti, Mende, Somali, Yoruba	Mallick et al. (2016)
Australasia	Australian, Bougainville, Hawaiian, Maori, Papuan	
Caucasus	Abkhasian, Adygei, Armenian, Chechen, Georgian, Lezgin, North Ossetian	
Central America	Mayan, Mixe, Mixtec, Pima, Zapotec	
Central Asia, Siberia	Aleut, Altaian, Chukchi, Eskimo Chaplin, Eskimo Naukan, Eskimo Sireniki, Even, Itelman, Kyrgyz, Mansi, Tajik, Tlingit, Tubalar, Ulchi, Uygur, Yakut	
East Asia	Ami, Atayal, Burmese, Cambodian, Dai, Daur, Dusun, Han, Hezhen, Igorot, Japanese, Kinh, Korean, Lahu, Miao, Mongola, Naxi, Oroqen, She, Thai, Tu, Tujia, Xibo, Yi	
Europe	Albanian, Basque, Bergamo, Bulgarian, Crete, Czech, English, Estonian, Finnish, French, Greek, Hungarian, Icelandic, Norwegian, Orcadian, Polish, Russian, Saami, Sardinian, Spanish, Tuscan	
North Africa	Mozabite, Saharawi	
South America	Chane, Karitiana, Piapoco, Quechua, Surui	
South Asia	Balochi, Bengali, Brahmin, Brahui, Burusho, Hazara, Irula, Kalash, Kapu, Khonda_Dora, Kusunda, Madiga, Makrani, Mala, Pathan, Punjabi, Relli, Sindhi, Yadava	
West Asia	BedouinB, Druze, Iranian, Iraqi Jew, Jordanian, Palestinian, Samaritan, Turkish, Yemenite Jew	
Ancient human	Ust’-Ishim (45000 YBP Siberian)	
Archaic human	Altai Neanderthal, Denisovan	
Tibetan/Sherpa	Genome	East Asia	Tibetan (AAGC022053, AAGC022056, AAGC022062, AAGC022070, AAGC022079, AAGC022109, AAGC022067), Sherpa (AAGC022072, AAGC022073, AAGC022074, AAGC022075, AAGC022076)	Lu et al. (2016)
Lazaridis et al. (2016)	Array	Africa	BantuKenya, Biaka, Esan, Gambian, Ju hoan North, Luhya, Luo, Mandenka, Masai, Mbuti, Mende, Somali, Yoruba	Lazaridis et al. (2016)
Australasia	Australian, Papuan	
Caucasus	Abkhasian, Adygei, Armenian, Chechen, Georgian, Lezgin, North Ossetian	
Central America	Mayan, Mixe, Mixtec, Pima, Zapotec	
Central Asia, Siberia	Aleut, Altaian, Chukchi, Even, Kyrgyz, Mansi, Tajik, Tubalar, Ulchi, Uygur, Yakut	
East Asia	Ami, Atayal, Burmese, Cambodian, Dai, Daur, Dusun, Han, Hezhen, Japanese, Kinh, Korean, Lahu, Miao, Mongola, Naxi, Oroqen, She, Thai, Tu, Tujia, Xibo, Yi	
East Eurasian	Dolgan, Eskimo, Itelmen, Kalmyk, Koryak, Nganasan, Selkup, Tuvinian, Yukagir	
Europe	Albanian, Basque, Bulgarian, Czech, English, Estonian, Finnish, French, Greek, Hungarian, Icelandic, Norwegian, Orcadian, Russian, Sardinian, Spanish	
North Africa	Mozabite, Saharawi	
South America	Karitiana, Piapoco, Quechua, Surui	
South Asia	Balochi, Bengali, Brahui, Burusho, Hazara, Kalash, Kusunda, Makrani, Pathan, Punjabi, Sindhi	
West Asia	BedouinB, Druze, Iranian, Jordanian, Palestinian, Turkish	
Ancient human	Mal’ta MA1, Afontva Gora-2 (AG2), Ust’-Ishim	
Skoglund et al. (2016)	Array	Southeast Asia	Malay, Vietnamese	Skoglund et al. (2016)
Island Southeast Asia	Ilocano, Kankanaey, Murut, Tagalog, Visayan	
Ancient Lapita	Lapita Tonga, Lapita Vanuatu	
The HUGO Pan-Asian SNP Consortium	Array	China	CN-CC (Zhuang), CN-GA (Han), CN-HM (Hmong), CN-JI (Jiamao), CN-JN (Jinuo), CN-SH (Han), CN-UG (Uyghur), CN-WA (Wa)	HUGO Pan-Asian SNP Consortium et al. (2009)
India	IN-DR (Proto-Austroloids), IN-EL (Caucasoids (may have admixture with Mongoloids)), IN-IL (Caucasoids), IN-NI (Mongoloid features), IN-NL (Caucasoids), IN-SP (Caucasoids), IN-TB (Mon- goloid features), IN-WI (Caucasoids), IN-WL (Caucasoids)	
Indonesia	ID-AL (Alorese), ID-DY (Dayak), ID-JA (Javanese), ID-JV (Javanese), ID-KR (Batak Karo), ID-LA (Lamaholot), ID-LE (Lembata), ID-ML (Malay), ID-MT (Mentawai), ID-RA (Mangga- rai), ID-SB (Kambera), ID-SO (Manggarai), ID-SU (Sunda), ID-TB (Batak Toba), ID-TR (Toraja)	
Japan	JP-ML (Japanese), JP-RK (Ryukyuan)	
Korea	KR-KR (Koreans)	
Malaysia	MY-TM (Proto-Malay), MY-BD (Bidayuh), MY-KS (Negrito), MY-JH (Negrito), MY-KN (Malay), MY-MN (Malay)	
Pacific	AX-ME (Melanesians)	
Philippines	PI-AT (Ati), PI-AG (Agta), PI-MA (Minanubu), PI-MW (Mamanwa), PI-AE (Ayta), PI-IR (Iraya), PI-UI (Filipino), PI-UN (Filipino), PI-UB (Filipino)	
Singapore	SG-CH (Chinese), SG-ML (Malay), SG-ID (Indian)	
Taiwan	AX-AM (Ami), AX-AT (Atayal)	
Thailand	TH-HM (Hmong (Miao)), TH-KA (Karen), TH-LW (Lawa), TH-MA (Mlabri), TH-MO (Mon), TH-PL (Paluang), TH-PP (Plang), TH-TK (Tai Khuen), TH-TL (Tai Lue), TH-TN (H’tin), TH-TU (Tai Yuan), TH-TY (Tai Yong), TH-YA (Yao) TW-HA (Chinese), TW-HB (Chinese)	
Japanese Archipelago Human Population Genetics Consortium	Array	Japan	Ainu, Ryukyuan, Japanese	Japanese Archipelago Human Population Genetics Consortium (2012)

We merged all datasets using the PLINK software (Purcell et al., 2007) and prepared five different merged datasets (Supplementary Table 2). F23 genotypes were then merged to each of the combined datasets. We downloaded the bam files of Mal’ta MA1, a 24000 years before present (YBP) Siberian genome (Raghavan et al., 2014) (http://www.cbs.dtu.dk/public/malta/MA1_AG2/bams/); Tianyuan, a 40000 YBP East Asian genome (Yang et al., 2017) (ENA: PRJEB20217); Devil’s Gate 1 and 2, 7700 YBP Neolithic East Asians from the Amur Basin (Siska et al., 2017) (ENA PRJEB14817); and Ikawazu Jomon IK002 and prehistoric Southeast Asians (McColl et al., 2018) (ENA PRJEB26721). We also downloaded the sra file of Chokhopani C1, a 3150–2400 YBP Nepal genome from the Annapurna Conservation Area (ACA) (Jeong et al., 2016b) (ENA SRP065070 and sample accession no. SRR3728363) and fastq files of Kostenki 14, a 38700–36200 YBP West Eurasian genome (Seguin-Orlando et al., 2014), and mapped them to hs37d5. For these genomic data and genomic data of our three Jomon individuals, F5, Sanganji131421-3, and Sanganji 131464, we called an allele by sampling a random read covering each SNP site with pileupCaller (https://github.com/stschiff/sequenceTools) used by Schuenemann et al. (2017), and merged these with the already merged datasets.

Admin
Administrator

Posts: 81,722

Jomon: Paleolithic Contingent in Modern Japanese Sept 22, 2021 22:49:19 GMT

Quote

Post by Admin on Sept 22, 2021 22:49:19 GMT

Principal component analysis (PCA)
SNPs observed in modern human populations were used to examine the genetic relationships between individuals through PCA. For genomic data, two filters, i.e., minor allele frequency <0.005 and genotype call rate <0.9, were applied using PLINK. Finally, we performed PCA using the smartpca program in the EIGENSOFT package (Patterson et al., 2006). When we included multiple ancient samples in the PCA, we used lsqproject option YES to project all the ancient samples. We used ancient samples having a genomic coverage of more than 2000 SNP coverage sites.

Phylogenetic analysis
Phylogenetic relationships between ancient and modern humans were first estimated using the neighbor-joining method (Saitou and Nei, 1987) based on Nei’s (1972) minimum distances. Trees were produced using SplitsTree4 software (Huson and Bryant, 2006).

For phylogenetic reconstruction with gene flow events, TreeMix (Pickrell and Pritchard, 2012) was used. Using Mbuti pygmies/Denisova as the root, trees were constructed assuming no migration and an incremental number of assumed gene flow events. Each run was performed in blocks of 500 SNPs to account for linkage disequilibrium, with sample size correction turned off.

F-statistics to detect gene flow
The f3-statistic (qp3Pop), the f4-statistic (qpDstat with f4 mode: YES), and the f4 ratio test, which are part of Admix-Tools (Patterson et al., 2012), can infer gene flows without invoking demographic models. These statistics are based on the number of shared alleles between test populations. We used |Z| > 3 as a significance cutoff.

We used the outgroup and admixture f3-statistics. The outgroup f3-statistic, f3 (Outgroup; X, Y), measures the shared genetic drift between two source populations, X and Y, compared with the outgroup population. Larger values indicate greater shared genetic drift between populations X and Y, which suggests a closer relationship. The admixture f3-statistic, f3 (Test; X, Y), detects whether the test population is explained by an admixture of populations X and Y. A significantly negative statistic implies that populations X and Y are candidates of the genetic source for the test population.

In the f4-statistic, f4 (O, X; Y, Z) measures allele frequency differences between two pairs of the group, (O, X) and (Y, Z), assuming population O as an outgroup to the three other populations. If the topology (O, X) (Y, Z) is correct, the frequency differences between Y and Z should reflect genetic drift, which is not correlated with those between O and X. Thus, the expected value of the product of frequency differences, f4, is zero. If f4 is significantly positive, X and Z share more alleles than Y and Z share, and if the f4 is significantly negative, X and Y share more alleles than X and Z share. To minimize the effects of possible biases due to small number of SNPs, only datasets with more than 50000 SNPs were used for this test.

To estimate the proportion of the Jomon ancestry in modern humans, we used the f4 ratio test, assuming the tree topology. There are at least three populations in the Japanese archipelago, the Ainu, Mainland Japanese, and Ryukyuan, who live in the northernmost, middle, and southernmost parts of the archipelago, respectively (Supplementary Figure 1). To reduce the influence of the recent admixture from Mainland Japanese into Ainu, only seven individuals who were most distant from the Mainland Japanese in the PCA plot were used in this analysis.

Testing models of phylogenetic relationships
To confirm the phylogenetic relationships and admixture events among populations illustrated by PCA, TreeMix, and F-statistics, we constructed the admix-tree using the qp-Graph program of the AdmixTools software. We chose five SGDP populations (Mbuti, Sardinian, Mixe, Han, and Japanese), MA1, and F23 for the graph construction, and we selected Mbuti, Sardinian, and Han as non-admixed populations. We then added F23, MA1, Mixe, and Japanese, sequentially. We calculated all possible admixed and non-admixed trees when we added the new population to the tree and kept the trees that produced the graph with the better Z-score (|Z| > 3) and fitting scores.

Estimating population size changes and divergence times
Changes in effective population size (Ne) over time and divergence times were estimated using the Multiple Sequentially Markovian Coalescent (MSMC) program (Schiffels and Durbin, 2014). Genotype calls by HaplotypeCaller described in the ‘Haplotype-based analysis’ section were used for MSMC analysis. For this analysis, we retained singleton variants but filtered out genotypes where genotype quality was less than 20 or where the fraction of sequence reads supporting the minor allele was less than 0.2. Haplotype phasing was conducted with Eagle version 2.3.2. Variants with depth between one-half and twice the average depth of each BAM file were used for the MSMC run. For estimation of Ne over time, four haplotypes (two individuals) representing each population were used, with the exception of F23 (one individual). Divergence times between populations were estimated in a pairwise manner, with one individual used to represent their respective population. For both of these analyses, a 30 year generation time and mutation rate of 1.25 × 10−8/site/generation were assumed. The mid-point of relative cross-coalescence rate (RCCR 0.5) was used for estimating the divergence time between populations.

Admin
Administrator

Posts: 81,722

Jomon: Paleolithic Contingent in Modern Japanese Sept 23, 2021 6:10:51 GMT

Quote

Post by Admin on Sept 23, 2021 6:10:51 GMT

Results
Mapping and DNA sequence authenticity
Both of our Funadomari Jomon samples showed a higher proportion of reads assigned to humans (77% and 50% in F23 and F5, respectively; see Figure 3), which are much higher than those found in most of other studies (e.g. Kanzawa-Kiriyama et al., 2017). Approximately 1.32 billion and 0.13 billion reads from DNA libraries for F23 and F5 were uniquely mapped to the human reference genome (GRCh37), and the peaks of depth were 48- and 1-fold, respectively (Figure 4a). We deposited all the short-read sequences produced for this study into the DDBJ Sequence Read Archive (SRA; trace.ddbj.nig.ac.jp/dra/). Sequence data are available for download through SRA (accession number DRA008031). Mapped data are available upon request.

Figure 3
Content and distribution of the metagenomes for F5 and F23. About 80% and 50% of the metagenomes were assigned to humans or primates with BLASTN, whereas bacterial DNA accounted for approximately 10% and 30% of sequences, respectively. About 10% and 20% of the reads were not classified into specific taxa.

Figure 4
Summary of mapped reads. (a) Distribution of read depths for F5 and F23. (b) Fragment length distribution after mapping to the reference genome, filtering with mapq30, collapsing PCR duplicates, removing clipped-reads, realigning, recalibrating, and filtering with read length.

The read length distribution for each library is shown in Figure 4b. The peak insert length of mapped sequence reads for F5 and F23 was 101 bp in both individuals, which is much longer than that of Sanganji Jomons, which is around 80 bp (Kanzawa-Kiriyama et al., 2017), indicating good DNA preservation. From the ratio of sequence reads mapped to the X and Y chromosomes (Table 2), we concluded that F5 and F23 were male and female, respectively. This identification is consistent with the morphological observations of these individuals (Matsumura and Anezaki, 2000). We observed characteristic features of ancient DNA degradation, including C to T misincorporation, at high frequencies in the 3′ and 5′ ends, and an increased fraction of purines close to the sequence read termini (Supplementary Figures 2 and 3), indicating that the sequence reads most likely contained endogenous Jomon DNAs. The ratio of the misincorporation was lower than those of previous studies, indicating that the DNA preservation was good compared to that of other ancient samples, probably because of the cold environment in the Funadomari site.