Genomic detection of a secondary family burial in a single jar coffin in early Medieval Korea

Abstract Objectives Family relationship is a key to understand the structure of past societies but its archeological reconstruction mostly stays circumstantial. Archaeogenetic information, especially genome‐wide data, provide an objective approach to accurately reconstruct the familial relationship of ancient individuals, thus allowing a robust test of an archaeology‐driven hypothesis of kinship. In this study, we applied this approach to disentangle the genetic relationship of early Medieval individuals from Korea, who were secondarily co‐buried in a single jar coffin. Materials and Methods We obtained genome‐wide data of six early Medieval Korean individuals from a jar coffin. We inferred the genetic relatedness between these individuals and characterized their genetic profiles using well‐established population genetics methods. Results Congruent with the unusual pattern of multiple individuals in a single jar coffin, genome‐wide analysis of these individuals shows that they form an extended family, including a couple, their two children and both paternal and maternal relatives. We show that these early Medieval Koreans have a genetic profile similar to present‐day Koreans. Discussion We show that an unusual case of the secondary multiple burial in a single jar coffin reflects family relationship among the co‐buried individuals. We find both paternal and maternal relatives co‐buried with the nuclear family, which may suggest a family structure with limited gender bias. We find the genetic profile of early Medieval Koreans similar to that of present‐day Koreans, showing that the genetic root of the present‐day Koreans goes back at least 1500 years in the Korean peninsula.


| INTRODUCTION
The burial practice using jar coffins was widespread throughout past societies in East Asia including Southeast Asia, northern China, Manchuria, the Japanese archipelago, and the Korean Peninsula, starting from the Neolithic period and continuing to the historic times (Bacvarov, 2006;Boeyens et al., 2009;W. Kim, 1973;Shewan et al., 2020). In the Korean peninsula, while prehistoric and historic burials with jar coffins are frequently found in all regions, they are most concentrated in the southwestern Korea, such as Naju and Yeongam regions, during the 4-6th centuries AD (E.-J. Kim, 2021; E. K. Kim, 2020;Kim et al., 2010;Oh, 2008;Park, 2010). Although burials with large, human-height-scale jar coffins are famous, most jar coffins found in Korea are much smaller than human body. Mostly, these small jar coffins host a single individual's skeleton secondarily collected after cremation. The practice of cremation and the following secondary burial often result in a loss of some skeletal elements, destruction of detailed morphological features that convey information on the sex, age, and pathological history of the individual, and poor preservation of biological macromolecules such as proteins and nucleic acids. For this reason, skeletal remains from the jar coffins have been investigated in bioanthropological and archaeogenetic studies only in a limited scale, leaving the characteristics of people buried in the jar-coffins an open question (M. J. Kim et al., 2010; K.-S. Lee et al., 1999).
The Dangbuk-ri archeological site is located at Gunsan city in the west coastal region of South Korea (Figure 1). Excavated during Summer 2016, this site includes a total of 21 burials of the early Medieval period (5-7th century AD; or the "three kingdoms" period). Among these 21 burials, 16 stone-cist and 4 stone-chamber burials did not yield skeletal remains, while the remaining one jar coffin burial housed skeletal remains of multiple individuals (Figure 2). Considering the burial context, it is most likely that the skeletal elements of these individuals were collected secondarily into this jar coffin and buried at the same time. This case is considered unusual because most jar coffin burials house only one individual and because individuals in the other types of multiple burials in the early Medieval Korea were usually added into the burial in a sequential manner over an extended period of time (E.-J. Kim, 2021). For example, in stone-chamber tombs with multiple individuals interred in a sequential manner, skeletal elements of individuals buried earlier tend to be out of their anatomical positions and placed at the periphery of the chamber because they were moved away from the center when the next individual was interred later (J. Lee et al., 2006).
In this study, we investigated the genetic relatedness between these co-buried individuals and their genetic profiles using genome-wide data.
In addition to revealing the familial relationship among them, we model their genetic profiles in terms of their relationship with ancient and present-day East Asian populations, especially with present-day Koreans and ancient populations near the Korean peninsula, in high-resolution by utilizing rich information in genome-wide data.

| Archeological context of the Gunsan jar coffin
The Dangbuk-ri site is located in Gunsan, Jeollabuk-do province on the Korean peninsula ( Figure 1). The site was salvage excavated for the construction of a railroad line in 2016 by the Central Institute of Cultural Heritage. The archeological investigation was conducted during the period of approximately 60 days, starting from July 2016. Geographically the site was on the southern hillside of Mt. Dottae, and it covered the top of a hill, the height of which was about 60 meters above sea level. A total of 21 Three Kingdoms period burials were found at the Dangbuk-ri site. Among those burials, 16 were the stone-cist type, four burials were the stone-chamber type, and the remaining one was a jar coffin. The stone-cist and stone-chamber tombs around the jar coffin burial were formed in the Three Kingdoms period given their burial structure and construction techniques. Moreover, typological analysis on the jar coffin suggests the technique of jar manufacture closely resembles that of the grave goods pottery from stone-chamber and stone-cist tombs. Specifically, the shape of the ceramic bowl out of stone-cist and stone-chamber tombs closely resemble the pottery made in China around the 6th century CE. Considering these archeological contexts, the jar coffin from the Dangbuk-ri site was contextually dated to the mid-6th century CE.
The jar coffin burial was located on the southeastern hillside from the top, surrounded by stone-cist and stone-chamber burials without any specific pattern. The jar was laid in a circular pit and two stones were placed beside the jar to fix the jar body. At the time of discovery, the jar was found slanted slightly. The mouth area of the jar has been broken without a cover of the coffin. In the jar coffin, multiple skeletons were found together ( Figure 2). No specific pattern was detected in the placement of skeletal elements within the jar coffin, although more skulls were placed in the upper part of the jar. No case of a sequential multiple burial (i.e. individuals were sequentially put to the jar coffin over an extended time period) has been reported for this type of jar coffins, although there are cases with those consisting of two jars. Also, individuals in this jar coffin were not separated into layers, a pattern expected for a sequential multiple burial. Thus, multiple skeletons in the jar coffin are considered as interred together approximately at the same time as a secondary burial. After finishing the excavation, a pottery as a grave good was found at the bottom of a circular pit, implying that ritual behavior was performed before laying the jar-coffin in state. The height of the jar, from the lower part of the jar to the top is 72.3 cm, and its diameter is 32.1 cm.

| Permission for destructive analyses of skeletal elements
The excavator of the Dangbuk-ri archeological site, Central Institute of Cultural Heritage, gave a permission to research the skeletal F I G U R E 2 The site map of the Dangbuk-ri site. The map shows the location of the jar coffin and the other 20 three kingdoms period tombs. Among the 16 stone-chamber and 4 stone-cist burials, four (two stone-chamber and two stone-cist) were located about 2 km away from the location on the map. The inset in the middle shows a wide view of the Gunsan jar-coffin in situ and the skeletal remains of the Gunsan jar-coffin individuals.
remains for the purposes of their scientific research. Specifically, the authors were given a permission for destructive sampling for the purpose of DNA sequencing, mass spectrometry of peptides, radiocarbon dating, and stable isotope analysis. The archeological investigation of the Dangbuk-ri site at Gunsan, South Korea was conducted by the Central Institute of Cultural Heritage, from July to August of 2016.
Prior to the full-scale excavation, a public announcement was made based on the Funeral Services Related Act (Section 2) in order to find possible descendants of the individuals found at the site. As no skeletal remains had been claimed 60 days after the announcement, the excavation was allowed to proceed, and with permission of the Central Institute of Cultural Heritage, the skeletal remains were moved to the Bio-anthropology Laboratory of Seoul National University for bioanthropological examination and ancient DNA analysis.

| Skeletal analysis
All skeletons were macroscopically examined by the two authors (E.J.W and C.L.J). In order to estimate the minimum number of individuals (MNI) in a jar, fragments were first identified by element type, and conjoining process was conducted to refit fragments from the same bone.
Identified elements were classified as either adult or subadult according to the degree of dental development, long-bone epiphyseal fusion, diaphyseal length, cranial size and cortical thickness. Bones with fully fused epiphyses and within the adult-sized range elements were classified as adults. Meanwhile, skeletal elements with unfused epiphyses and within the size range of the subadults were classified into the subadult category. After the skeletal elements were sorted to individuals, age estimation was further refined based on morphological features.
The value of MNI was derived by sorting elements into lefts and rights, and then taking the greatest number as the final estimate following a published protocol (T. E. White, 1953). Severely fragmented skeletal elements such as ribs, vertebrae and phalanges were excluded from calculating MNI due to uncertainty in siding. Then, visual pair-matching was conducted to decide if principal limb bones were from a single individual through the comparison of right and left elements. Finally, different parts of the skeletal elements were segregated into individuals based on the examination of degenerative changes, articulation, robusticity, age, and sex, along with osteometric sorting (Adams & Byrd, 2014).
The sex and age were estimated for each skeletal individual, based on the Standards for Data Collection (Buikstra & Ubelaker, 1994;Ubelaker, 1999). Sex was estimated using morphological features of the skull and robusticity of limb bones. We were unable to use the pelvis for sex estimation because most pelvic elements in the jar-coffin were severely fragmented. The age at death of each individual was estimated by dental attrition, number of antemortem tooth loss, cranial suture closure, degeneration of the auricular surface and pubic symphyseal surface, and degenerative changes of joint in limb bones and the vertebral column. The age of subadults was estimated using the degrees of tooth formation and eruption, epiphyseal fusion or diaphyseal length (Buikstra & Ubelaker, 1994; T. D. White & Folkens, 2005).

| Sampling of skeletal materials
Sample selection was performed by E.J.W. and C.L.J. in the Bioanthropology laboratory (Department of Anthropology, Seoul National University). As compact bones are preferred for ancient DNA analysis, intact petrous parts of the temporal bone from seven individuals were chosen. Two out of nine individuals were not included in the analysis since the temporal parts including petrous were not well preserved. In case of five adults, the left petrous bone was selected, for two subadults, the right petrous parts were sampled. The outer surface of the petrous pyramid inside the skull was ground with a 4.8 mm cutter burr attached to a Dremel 9100-21 Fortiflex 2.5-Amp Stationary Flex Shaft Precision Rotary Tool (Dremel, Mount Prospect, IL). Then, the inferior border of the cochlea was ground to create a small opening into the osseous labyrinth (Sirak et al., 2017). Into this opening, a 3.2 mm engraving cutting burr was applied in a circular motion to obtain bone powder. The powder was collected in a sterilized paper foil and placed in a sterile 1.5 ml Eppendorf tube for DNA extraction.

| Ancient DNA laboratory work and sequencing
For each of the seven individuals, a double-strand double-indexed Illumina sequencing library was built from metagenomic DNA extracted from 30 to 50 mg of bone powder. DNA extraction and library preparation were performed using the previously published protocols (Dabney et al., 2013). For library preparation, a partial treatment of the uracil-DNA-glycosylase (UDG) enzyme was included following a published protocol  to confine deaminated bases to the ends of the reads. All laboratory works up to library prepraration step were performed in a dedicated ancient DNA (aDNA) clean room facility of the Max Planck Institute for the Science of Human History (MPI-SHH), Jena, Germany. For six of the seven individuals with sufficient levels of human DNA preservation, ranging 0.1%-2.6%, two rounds of insolution capture for 1.24 million ancestry-informative single nucleotide polymorphisms (SNPs) (the "1240K" panel) was performed to enrich libraries (Mathieson et al., 2015). All libraries were sequenced using single-end 76 base pair (bp) sequencing on the Illumina HiSeq 4000 platform following the manufacturer's protocols.
2.6 | Ancient DNA sequencing data processing and authentication aDNA sequencing data were processed using the EAGER v1.92.50 wrapper (Peltzer et al., 2016). Within the EAGER wrapper, Illumina adapter sequences were first trimmed from raw reads using the Adap-terRemoval v2.3.0 (Schubert et al., 2016). Adapter-trimmed reads of 30 bp or longer were mapped to the human reference genome with decoy sequences (hs37d5) using the aln/samse modules in BWA v0.7.17 with "-n 0.01" option (Li & Durbin, 2009). PCR duplicates were removed using the DeDup v0.12.5 (Peltzer et al., 2016), assuming that both ends of the reads were known. Unique mapped reads with Phred-scaled mapping quality score 30 or higher were kept using the samtools v1.9 . To remove deaminationbased misincorporations (5 0 C > T and 3 0 G > A), the first and last two bases of each read were soft-masked using the trimBam module of bamUtils v1.0.14 (Jun et al., 2015). Finally, pseudo-haploid genotype data were determined by randomly sampling a single high-quality base (Phred-scaled base quality score 30 or higher) per site per individual using samtools mpileup and pileupCaller v1.4.0.5 (downloaded from https://github.com/stschiff/sequenceTools). For C/T and G/A SNPs, end-masked BAM files were used. For the remaining SNPs that are not affected by post-mortem deamination, BAM files without endmasking were used to maximize sequence data usage. The authenticity of sequence data was checked by multiple measures. First, chemical modifications typical of aDNA molecules were tabulated using mapDamage v2.0.9 (Jonsson et al., 2013). Second, mitochondrial DNA contamination was estimated using the schmutzi v1.5.5.5 (Renaud et al., 2015). Third, X chromosome-based estimation of nuclear DNA contamination was performed for four male individuals using the contamination module of the ANGSD v0.929 program . Mitochondrial haplogroups were determined by applying the HaploGrep v2 program (Weissensteiner et al., 2016) to the consensus sequences called by the log2fasta program in schmutzi with "-q10 00 filter. Y haplogroups were determined using a modified version of the yHaplo program with "--ancStopThresh 10" option to prevent the root-to-tip search from halting in an internal branch due to missing data (downloaded from https://github.com/alexhbnr/yhaplo) (Poznik, 2016).

| Reprocessing of whole genome sequences of present-day Koreans
We downloaded high-coverage whole genome sequencing data of 104 Koreans from the Ulsan city from KoVariome data (ftp:// biodisk.org/Release/KPGP/KPGP_Data_2017_Release_Candidate/).
We kept properly aligned paired-end reads by applying "-f 0Â0003" filter in samtools view v1.9 ) and merged per-lane BAM files into per-individual ones using samtools merge. We removed duplicates using the Picard MarkDuplicates v2.20.0 (downloaded from https://broadinstitute.github.io/picard/) and kept reads with Phred-scaled mapping quality score 25 or higher. From these analysis-ready BAM files, we produced two genotype calls for the 1,233,013 SNPs in the 1240 K panel. First, we calculated genotype likelihoods for each individual using the UnifiedGenotyper module of the Genome Analysis ToolKit (GATK) v3.8.1.0 (McKenna et al., 2010) with the "--allSitesPLs" flag, calculated posterior genotype probability by multiplying genotype likelihoods with the GATK default prior [0.9985, 0.0010, 0.0005], and took genotype calls with posterior probability 0.900 or higher. Second, we randomly sampled a base with the Phred-scaled base quality score 30 or higher from a read with mapping quality score 30 or higher per site per individual to create a pseudo-haploid call that mimics the genotype calling procedure of low-coverage ancient individuals. For this procedure, we used samtools mpileup and pileupCaller v1.4.0.5. We merged the two genotype calls of present-day Ulsan Koreans with the genomewide data of world-wide present-day populations typed on the Affymetrix Axiom ® Genome-Wide Human Origins 1 array ("HumanOrigins") and the 1240 K dataset for downstream analyses. We determined genetic sex of each individual by comparing coverage of sex chromosomes with autosomes calculated for the 1240K sites using samtools depth ( Figure S1). We detected genetic outliers by projecting Ulsan Koreans to the top principal components calculated for 2077 present-day Eurasian individuals using the smartpca program v16000 in the EIGENSOFT package v7.2.1 (Patterson et al., 2006). For the remaining individuals, we detected close relative pairs by calculating pairwise mismatch rate (PMR) for each pair using the random haploid calls and by estimating kinship coefficients using the "--Z-genome" module in PLINK v1.90b6.9 (Chang et al., 2015). We removed one individual from each relative pair up to the second-degree relatives for the downstream group-based analyses ( Figure S2 and Table S1).
We provide a list of analysis groups and individuals used for each analysis (Table S2).

| Genetic kinship analysis
We calculated PMR between each pair of ancient Gunsan jar coffin individuals across the 1,150,639 autosomal SNPs in the 1240K panel.
Each pair was covered by at least 66,442 SNPs. We estimated probability of sharing 0, 1, and 2 alleles using the lcMLkin v0.5.0 to distinguish between parent-offspring and full sibling (Lipatov et al., 2015).
For the group-based analyses, we removed first-degree relatives and kept three individuals (GUC002, GUC003, and GUC005) with minimal genetic relatedness (Table S3).

| F-statistics and qpWave/qpAdm analysis
We obtained f-statistics using 1240K dataset to maximize SNP coverage of ancient individuals. We calculated outgroup-f 3 statistics of the form f 3 (Mbuti; Gunsan jar coffin, world-wide) using qp3Pop v650 to measure genetic affinity between two populations. We calculated f 4 statistics of the form f 4 (Mbuti, world-wide; Gunsan jar coffin, presentday Korean) using the qpDstat v970 with the option "f4mode: YES." We used 121 present-day world-wide populations and 167 ancient populations for these calculations (Table S2). We tested various admixture models of ancient and present-day East Asian populations using qpWave v1200 and qpAdm v1201 programs in the AdmixTools v7.0 (Lazaridis et al., 2016;Reich et al., 2012) on 1240K dataset. We  south genetic cline. And finally Funadomari Jomon is to tag the Jomon ancestry. We included "allsnps: NO" option.

| Ancient genome-wide data production
We performed a genetic investigation of early Medieval individuals from Dangbuk-ri, Gunsan, in the southwestern region of the Republic of Korea ( Figure 1 and Table 1). These individuals are dated to 5th to 7th century AD based on archeological contexts and are found in a single jar coffin, with signatures of secondary burial: multiple burials, either primary or secondary, within a single jar coffin are exceptional given that most jar coffins host only a single individual. Based on the macroscopic examination of all skeletal elements retrieved from the jar coffin, we determined MNI as nine, including six adults and three subadults (Table 1). Following an in-solution capture protocol for~1.2 million informative SNPs (Mathieson et al., 2015), we obtained genome-wide data with 189K-764K on-target SNPs covered at least once for six individuals ( Table 2) (Table 2). For each individual, we called haploid genotypes across the 1240K panel SNPs by randomly sampling a high-quality base. We concatenated their genotype data with published presenet-day and ancient individuals on the Huma-nOrigins and the 1240K data sets for downstream analyses (Table S2).
T A B L E 2 A summary of sequencing and genetic information of six ancient individuals in this study

| A familial relationship of individuals from a single jar coffin
We first ran PCA of 2077 present-day Eurasians and projected ancient individuals from the Gunsan jar coffin onto the top PCs ( Figure 3). All six individuals fall close to each other and to presentday Koreans, suggesting no substantial genetic heterogeneity among them and an overall close relationship with present-day Koreans ( Figure 3). In the PCA of present-day East Asians including presentday Koreans, the Gunsan jar coffin individuals also overlap with present-day Koreans ( Figure S3).
Then we measured the genetic relatedness between these ancient individuals to test if the unusual occasion of multiple individuals in a single jar coffin reflects their close relationship. Our estimates of genetic relatedness based on genome-wide data indeed show that these individuals form a closely related extended family: among 15 pairs combined from six individuals, we observe six first-degree, two second-degree, and one third-degree relative pairs ( Figure 4a; Table S3). Incorporating mitochondrial and Y haplogroup information and age at death, we distinguish between parent-offspring and full sibling pairs and propose the most plausible pedigrees (Figure 4b). At the core of this pedigree is a quartet family composed of a couple and their two children, one adult male and one subadult female. Among the remaining two individuals, one male (GUC005) is the first degree relative of the mother of the quartet (GUC004) and the second degree relative of the two children (GUC001 and GUC007), but is unrelated to the father (GUC002). The 1st degree relationship of GUC005 and GUC004 is likely a parent-offspring one given the near-zero sharing of both alleles at the same SNP (Table S3). Both father-daughter and mother-son relationship are compatible with genetic data, however, given that GUC005 and the quartet children share the identical mitochondrial haplogroup, we propose mother-son may be more likely.
That is, GUC005 and (GUC001, GUC007) may be half-siblings with the same mother. Lastly, GUC003 is most likely a 3rd degree relative of the quartet father, sharing the same Y haplogroup.
We also investigated the distribution of ROH segments in the Gunsan individuals using hapROH (Ringbauer et al., 2021). Finding few long ROH segments (>4 cM), we conclude that none of these individuals were the offspring of close relatives ( Figure S6).

| The genetic profile of early medieval Koreans
To understand the genetic profile of the early Medieval Koreans from the Gunsan jar coffin, we first measured the genetic affinity between the Gunsan jar coffin group and world-wide present-day and ancient populations using outgroup-f 3 statistics of the form f 3 (Mbuti; Gunsan jar coffin, world-wide). As expected, the Gunsan jar coffin group shows the highest genetic affinity with present-day Koreans ( Figure S4). Formally testing the genetic symmetry of the early Medieval and present-day Koreans using f 4 statistics of the form f 4 (Mbuti, world-wide; Gunsan jar coffin, present-day Korean), only a few groups break the symmetry without clear geographic distribution ( Figure S5).  (Table S6) Table S7A). In both ancient and present-day Koreans, we do not detect a statistically significant contribution from the Jomon hunter-gatherer gene pool of the Japanese archipelago (Table S7A), although previous studies report occasional presence of the Jomon ancestry contribution from Neolithic to the early Medieval period (Gelabert et al., 2022;Robbeets et al., 2021).
When we replace the genetic northern proxy from WLR_BA to Middle Neolithic individuals from the Miaogizou site in Inner Mongolia ("Miaozigou_MN"), we detect a small but significant amount of Jomon contribution in the Gunsan individuals and present-day Ulsan Koreans (3.1%-4.4%; Table S7B). We believe that WLR_BA provides a more suitable model for ancient and present-day Koreans given its geographical and temporal proximity to them. The remaining wellfitting source pairs provide qualitatively similar results (Table S8).

| DISCUSSION
We present a genomic study of individuals from the Gunsan jar coffin, an unusual case of a secondary multiple burial where at least nine individuals were co-buried in a single small jar coffin. We confirm our hypothesis that this unusual burial represents an unusual relationship among the co-buried individuals, that is, an extended family, including a core of a couple and their two children, as well as both paternal and maternal relatives. The inferred pedigree may imply little gender bias in the family structure of early Medieval Koreans lived in the area. Further archaeogenomic studies on the ancient individuals excavated from an unusual burial context or those from a single cemetery will provide more insights into past mortuary practices and social structures.
While the genetic origins of present-day Koreans have not been fully understood due to lack of relevant ancient genome data, individuals from the Gunsan jar coffin provide among the first past human genetic profiles in ancient Korea. Since this type of small jar coffins are considered to be associated with commoners rather than with  (Table S7A). For Kofun that has a substantial Jomon ancestry, we show the results of a threeadmixture model of WLR_BA+Xitoucun+Jomon_Ikawazu. Horizontal bars represent the standard error measure (s.e.m.) of ancestry proportion estimates, calculated by the 5 cM block jackknifing procedure as implemented in the qpAdm program.