Original Articles
Mitochondrial DNA sequence analysis of ancient human remains found in a 2000-year-old elite Xiongnu cemetery in northeast Mongolia
2023 Volume 131 Issue 1 Pages 55-64
Details
2023 Volume 131 Issue 1 Pages 55-64
The Xiongnu empire was a dominating nomadic tribe in the east of southern Siberia between the third century BC and first century AD. Archeological sites and artifacts of the Xiongu have been found primarily in Mongolia and scattered in Russia and the northern part of China. The historical sites are mainly tombs. Rectangular terrace tombs with satellite graves—small circular tombs around the main large terrace tombs—have been found throughout Xiongnu cemeteries. This is the first study to report a genetic analysis of the relationship among individuals in these satellite tombs. The ancient DNA samples were extracted following a strict decontamination protocol in a dedicated laboratory using a silica-based method. The mitochondrial DNA segments required for haplogrouping were assessed using multiplex and nested real-time polymerase chain reaction-based methods to amplify multiple target segments simultaneously. This reduced the consumption of valuable samples, labor, and time. Mitochondrial DNA haplotypes were analyzed from the skeletons excavated from six satellite tombs found in a 2000-year-old Xiongnu elite cemetery at Duurlig nars, Mongolia. Interestingly, the different mitochondrial haplotypes (A, C4, D4, G1, G2, W) of six satellite tombs did not reveal any maternal kinship among the six individuals in the tombs. Two sacrifice tombs had two females in their twenties within them. Future studies should assess Y haplotypes, Y chromosomal short tandem repeats, autosomal short tandem repeats, and phenotype and biogeography ancestry-informative DNA.
The Xiongnu empire ruled over a wide territory, from Inner Mongolia to Lake Baikal and from Liaoning to parts of Xinjiang (Psarras, 1995). Xiongnu historical sites generally consist of tombs, with some ruined castles, settlements, and workshop sites, including iron manufacture ruins or kiln sites (Oh et al., 2018). Xiongnu tombs are widespread throughout the Eurasian steppes, including Mongolia, Russia, China, and Kyrgyzstan. Xiongnu tombs were constructed at different periods according to historical changes and migration (Lim et al., 2016). The population of Xiongnu was likely composed of various ethnicities due to its vast territory (Kim et al., 2010).
Xiongnu tombs in Mongolia are divided into two different types: rectangular terrace and circular tombs. They can be identified by a stone line on the surface of the tombs. The rectangular terrace tombs are upside down pyramidal shapes with a dromos and the entry on the south side. This appearance gives an external convex shape. Most of the terrace tombs with a dromos are quite large and considered as Xiongnu elite tombs (Chang et al., 2007; Yang and Eregzen, 2017). The large terrace tombs of Xiongnu elites are found in nine regions of Mongolia: Takhiltyn khotgor, Ovoono khar, Gol mod I, Gol mod II, Khyalganat, Noyon uul, Selbe, Bor bulag, and Duurlig nars. It was hypothesized that the convex-shaped tombs of Xiongnu elites included the highest classes based on their size and the first-class artifacts found within them (Yang and Eregzen, 2017). Molecular genetic studies have been performed on ancient human skeletons in the tombs of different archeological sites to answer various anthropological questions (Kaestle and Horsburgh, 2002; Kakuda et al., 2016; Hollard et al., 2018). Furthermore, kinship between the ancient human skeletons has been analyzed in Xiongnu cemeteries (Keyser-Tracqui et al., 2003; Keyser et al., 2021). There are several opinions about the relationship between the Xiongnu terrace tombs and their satellite tombs; however, no study has elucidated the molecular genetics.
We analyzed the mitochondrial haplogroup of six ancient human skeletons excavated from satellite tombs in a 2000-year-old Xiongnu elite cemetery at Duurlig Nars, Mongolia to study the relationship between these tombs. To facilitate mitochondrial DNA (mtDNA) sequencing analysis of these ancient samples, we employed real-time polymerase chain reaction (PCR)-based multiplex and nested PCR methods.
The Xiongnu Elite cemetery, dating back to 2000 BP was discovered at Duurlig Nars (N 48° 32' 40.9", E 111° 04' 42.9") in the Bayan Adarga sum (district) north of Khenti Aimag, Mongolia (Figure 1). Khenti Aimag is the birthplace of Genghis Khan. It has been described previously (Kim et al., 2010). In brief, Duurlig Nars is 450 km northeast of the capital of Mongolia, Ulaanbaatar. Over 280 tombs were found in the Xiongnu cemetery (Figure 2A).
Location of Duurlig Nars of Mongolia (Google maps)
(A) Map of the Duurlig Nars Xiongnu elite cemetery. The tombs in the open circle were excavated. (B) The main Xiongnu elite tomb in the open circle in (A) is 55 m in length, with a 32 m dromos. The satellite tombs around the large elite rectangular tomb are located on the east side (E1, E2, E3), on the south (S1, S2), and on the west (W1, W2, W3, W4, W5, W6).
Korean–Mongol archeologists excavated four tombs, including three satellite tombs, in Duurlig Nars during 2006–2007. The molecular genetic results of three human skeletons from three satellite tombs were analyzed (Kim et al., 2010). During 2010–2011, the large square tomb complexes (a rectangular tomb with satellite tombs) were excavated. The large pillaged tomb was an Xiongnu elite tomb with a convex shape (55.5 m total length, 32 m dromos or entry chamber, and 15 m depth). In total, 11 satellite tombs were excavated: three in the east, six in the west, and two in the south (Figure 2B) (National Museum of Korea et al., 2014).
Ancient human samplesThe samples used in this study are presented in Table 1. The bone dating was based on previous archaeological evidence (Kim, 2014).
| No. | Code | Excavation site | Age | Sample |
|---|---|---|---|---|
| 1 | MNE1 | Duurlig Nars, Mongolia | 100 BC–100 AD, Xiongnu | Femur |
| 2 | MNE2 | Duurlig Nars, Mongolia | 100 BC–100 AD, Xiongnu | Femur |
| 3 | MNE3 | Duurlig Nars, Mongolia | 100 BC–100 AD, Xiongnu | Tibia |
| 4 | MNW1 | Duurlig Nars, Mongolia | 100 BC–100 AD, Xiongnu | Femur |
| 5 | MNW3 | Duurlig Nars, Mongolia | 100 BC–100 AD, Xiongnu | Humerus |
| 6 | MNW4 | Duurlig Nars, Mongolia | 100 BC–100 AD, Xiongnu | Molar teeth |
We strictly followed previously described precautions for ancient DNA (aDNA) extraction (Kim et al., 2008). All procedures were conducted in four separated and restricted clean rooms dedicated to aDNA experiments and under sterile conditions, using UV-irradiated gowns with head and leg covers, latex gloves, and face- and mouth-masks. The first room was for bone fragmentation and bone surface removal. The second was for bone surface decontamination, powdering, and DNA extraction. This room was maintained under HEPA-filtered positive air pressure to minimize potential contamination from the outside. The third was for PCR preparation and analysis. The fourth was for post-PCR analysis. We ensured that there was no sharing of materials and instruments between rooms. Two-hour-autoclaved or human-DNA-free certified materials were used. Workplaces, containers, and bench surfaces were cleaned by undiluted bleach and rinsing with DNA-free water and UV irradiation at 254 nm for at least one day. All steps were performed in a fume hood and laminar airflow clean bench. Extraction and template blanks were included through the entire process. Human-DNA-free aerosol-barrier pipettes and pipette tips were used throughout all the manipulations. For each sample, different bone parts or teeth were analyzed to ensure consistency in the results.
Ancient bone pretreatment and decontaminationWe followed our established protocol to eliminate potential exogenous DNA and soil-originated contaminants on the bone surface that can cause false results from potentially contaminated contemporary human DNA or analytic failure caused by severe PCR inhibition.
Bone fragmentation and surface removalA bone sample was cut at the thickest and hardest cortex area into several fragments using a cut-off wheel and an electric rotary tool (RTX-1, Black & Decker) in the fume hood. The entire outer and inner surface of the bone fragments was removed to a depth of at least 1 mm using sanding bands to remove humic substances, which are the primary inhibitory compounds of DNA extracted from soil or buried biological remains (Tebbe and Vahjen, 1993; Watson and Blackwell, 2000; Sutlovic et al., 2005, 2008). The bone fragments were further cut into rectangular pieces approximately 5 mm in size.
DecontaminationOne to two grams of the bone fragments were immersed in about 40 ml of undiluted bleach in a 50 ml centrifuge tube shaking intermittently for 30 min and rinsed in the tube, shaking gently with 500 ml of sterile DNA-free distilled water. The bleach treatment is effective for the destruction of contaminated human DNA on the bone surface (Kemp and Smith, 2005). The bone fragments were dried under UV irradiation of all surfaces overnight on a clean bench.
DNA extractionThe DNA extraction steps consisted of bone powdering, lysis of the bone powders, preparation of silica suspension in water, silica capture and washing of DNA from bone lysate, and DNA elution from silica particles. We followed a previously published method for DNA extraction from bones (Rohland and Hofreiter, 2007) with some modifications. We modified the washing and elution steps of DNA-bound silica particles, where we used buffers included in Qiaquick PCR purification kits and added a concentration step where the Qiaquick PCR purification kit was used.
Bone powderingThe decontaminated bone fragments were placed into a mixer mill jar and incubated in liquid nitrogen for 5 min. Bone powders were made with Mixer Mill MM301 (Retsch) for 1 min at full power. The fine bone powders were collected into a leakproof 50 ml centrifuge tube (Orange Scientific) and were temporarily stored in the refrigerator until used for DNA extraction.
Lysis of bone powdersThe bone powders (average 1.9 g; 1.6–2.5 g) were incubated with 10 ml lysis solution (0.45 M EDTA and 0.5 mg/ml proteinase K) at 30°C for 12 hours under rotation. The supernatant was collected in another 50-ml centrifuge tube after centrifugation at 3500 g for 10 min. The remaining pellets were saved in the freezer for a possible repeated lysis in the future (Table S1).
Preparation of silica suspension in waterSilica suspension was prepared for use in the capture of DNA in a lysate solution, as follows. First, 4.8 g of silica powder (silicon dioxide, Sigma S5631) was suspended in the distilled water to a final volume of 40 ml in a 50 ml centrifuge tube. After sedimentation for 1 hour, 39 ml were transferred into a new tube and left to settle for an additional 4 hours. Thereafter, 35 ml was removed and 48 μl of fresh concentrated HCl was added to the sediment to adjust the pH. After the silica suspension was divided into 400 μl aliquots in screwcap microtubes, the tubes were tightly closed, autoclaved for 2 hours, and stored in a refrigerator in the dark. The silica suspension was kept for use for up to 2 months.
Silica capture and washing of DNA from bone lysateThe bone lysate (~10 ml) was mixed with 30 ml binding buffer (5.5 M guanidine thiocyanate, 0.025 M sodium citrate, pH 7.0), 10 ml isopropanol, 90 μl of the silica suspension pre-prepared above, and 700 μl 10 M HCl in a 50 ml centrifuge tube. The mixture was incubated for 3 hours at room temperature under rotation. The DNA-bound silica particles were pelleted by centrifugation at 3500 g for 5 min. Most of the supernatant was discarded but some (<1 ml) was retained for the resuspension. The silica pellets were resuspended with the residual supernatant and transferred into a microtube. The silica particles were pelleted in the microtube by centrifugation at 5000 g for 1 min. The supernatant was removed, and the pellet was completely resuspended by pipetting up and down with 0.8 ml PE buffer (DNA wash buffer; QIAquick PCR purification kit, Qiagen). The PE-silica suspension was centrifuged and the supernatant was discarded. The silica pellet was resuspended in PE buffer. This wash step was repeated again. The silica pellet was collected and the residual PE buffer in the tube was completely removed from the pellet using a micropipette. The silica pellet in the microtube was dried with the lid slightly opened for 15 min in a dry bath at 50°C to completely evaporate the alcohol derived from the PE buffer.
DNA elution from silica particlesThe dried silica pellet was resuspended in 300 μl EB buffer (DNA elution buffer; QIAquick PCR purification kit, Qiagen) and incubated for 2 min at room temperature. The suspension containing the released DNA and silica particles was transferred into a 1.5 ml microtube and centrifuged, and only the supernatant was carefully collected in a new tube without disturbing the silica pellet. The supernatant (~250 μl) was concentrated using a QIAquick PCR purification kit (Qiagen) with a 50 μl elution. This elution was used as the aDNA template for subsequent DNA analysis such as PCR and sequencing.
Mitochondrial DNA analysis Primer designMultiplex primers were designed for the simultaneous amplification of two HV1 and two HV2 region segments in the control region. In addition, primers were designed for haplogroup-determining segments in the coding region. The most conserved primer positions for target segments were searched using the Haplotracker application (Kim et al., 2022) with the following inputs: range, an approximately 400-bp region including the target segments; haplogroup target(s), none for HV region segments and haplogroup for coding region segments. Acceptable primer cross-complementarities for multiple primers, melting temperature (Tm), and amplicon size were analyzed using LightCycler Probe Design 2.0 software (Roche). The primer Tms were designed to be less than 60°C for the control region and less than 65°C for the coding region. The amplicons were smaller than 300 bp. Primer specificity in Homo sapiens was checked using primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). For each target region, a pair of primers for the first round PCR and another pair for the second (nested PCR) were designed (Table 2).
| Segmenta | Primer nameb | Sequence (5'→3') | Tm (°C)c | Size (bp) |
|---|---|---|---|---|
| Control region | ||||
| 16003–16196 (HV1) | F15969 | CTTTAACTCCACCAYTAGCAC | 60.7–62.4 | 287 |
| R16255 | CTTTGGRGTTGCAGTTGATG | 61.8–64.2 | ||
| FN15986 | GCACCCAAAGCTAAGAT | 59.0 | 231 | |
| RN16216 | TTGCTGTACTTGCTTGTAAG | 59.7 | ||
| 16204–16413 (HV1) | F16186 | CCCYCCCCATGCTTACAA | 62.4–65.1 | 270 |
| R16455 | CGGAGCGAGGAGAGTA | 60.5 | ||
| FN16197 | CTTACAAGCAAGTACAGCAA | 59.4 | 236 | |
| RN16432 | TGTGCGGGATATTGATTTC | 58.9 | ||
| 35–263 (HV2) | F16565 | CGATGGATCACAGGTCTAT | 60.7 | 294 |
| R289 | TTTTGTTATGATGTCTGTGTGG | 60.4 | ||
| FN17 | CCCTATTAACCACTCACG | 59.1 | 265 | |
| RN281 | TGATGTCTGTGTGGAAAG | 58.2 | ||
| 184–392 (HV2) | F154 | TTATTTATCGCACCTACGTTCA | 60.0 | 268 |
| R421 | AGTGCATACCGCCAAAAG | 63.5 | ||
| FN163 | GCACCTACGTTCAATATTACA | 58.8 | 259 | |
| RN412 | CGCCAAAAGATAAAATTTGA | 56.2 | ||
| Coding regiond | ||||
| 5051–5167 (G) | F5007 | GCATACTCCTCAATTACCCACATAG | 64.2 | 213 |
| R5210 | GGTGGATGGAATTAAGGGTGTTAG | 64.5 | ||
| FN5026 | ACATAGGATGAATAATAGCAGTTCT | 61.9 | 160 | |
| RN5191 | GTTAGTCATGTTAGCTTGTTTCAG | 61.6 | ||
| 15280–15390 (G1) | F15242 | GGAGGCTACTCAGTAGACAGT | 64.3 | 196 |
| R15437 | CGAGGGCGTCTTTGATTGT | 65.2 | ||
| FN15264 | CCACCCTCACACGATT | 61.1 | 147 | |
| RN15410 | GGTGGAAGGTGATTTTATCG | 60.4 | ||
| 13507–13605 (G2) | F13463 | GCAGCCTAGCATTAGCAGGAATA | 65.1 | 186 |
| R13648 | GGAAGCGAGGTTGACCTGTTA | 64.9 | ||
| FN13486 | CCTTTCCTCACAGGTTTCTAC | 62.0 | 142 | |
| RN13627 | GGGTGAGAAGAATTATTCGAGT | 61.7 | ||
| 1690–1815 (A) | F1652 | CTTAACTTGACCGCTCTGAGCTAAAC | 66.9 | 201 |
| R1852 | GCAGAAGGTATAGGGGTTAGTCCTTG | 67.3 | ||
| FN1670 | AGCTAAACCTAGCCCCAAAC | 64.5 | 171 | |
| RN1840 | GGGGTTAGTCCTTGCTATATTATGC | 64.0 | ||
| 8408–8503 (D4) | F8358 | ACAGTGAAATGCCCCAACTAAATAC | 65.5 | 190 |
| R8547 | AGCGAACAGATTTTCGTTCATTTTG | 65.4 | ||
| FN8388 | TATGGCCCACCATAATTACC | 60.2 | 140 | |
| RN8527 | TTTTGGTTCTCAGGGTTTGTTATA | 60.6 | ||
| 15728–15901 (W) | F15699 | GCCCACTAAGCCAATCACTT | 64.7 | 235 |
| R15933 | CATCTCCGGTTTACAAGACTGG | 64.8 | ||
| FN15706 | AAGCCAATCACTTTATTGACTC | 61.1 | 221 | |
| RN15926 | GGTTTACAAGACTGGTGTATTAGTT | 62.4 | ||
| 3493–3598 (C) | F3453 | CGCTGACGCCATAAAACTCT | 65.4 | 188 |
| R3640 | CTAGGCTAGAGGTGGCTAGAAT | 64.7 | ||
| FN3473 | TCACCAAAGAGCCCCTAAAA | 61.5 | 146 | |
| RN3618 | AAATAGGAGGCCTAGGTTGA | 61.5 | ||
| 15133–15207 (C4) | F15076 | AGAAACCTGAAACATCGGCATTAT | 65.5 | 195 |
| R15270 | AGGGTGGGACTGTCTACTGA | 65.5 | ||
| FN15113 | ACTATAGCAACAGCCTTCAT | 61.7 | 115 | |
| RN15227 | CTAGGTCTGTCCCAATGTAT | 60.6 |
HV, hypervariable region; Tm, melting temperature.
a Numbers indicate nucleotide positions based on the revised Cambridge reference sequence (Andrews et al., 1999).
b F, forward primer; R, reverse primer; N, primers for nested PCR and sequencing.
c Primer melting temperature was calculated by using LC PDS software v. 2.0.
d Determining haplogroup is shown in parentheses.
To analyze mtDNA sequences of the aDNA extracts, the HV1 and HV2 segments of the control region were first subjected to PCR amplification, and then haplogroup-specific coding-region segments were targeted based on the sequencing results. A LightCycler 96 (Roche) system was used for the real-time multiplex PCR amplification of the four control region segments. We employed a nested PCR strategy (Kim, 2010), in which the first-round multiplex PCR targeting multiple mtDNA segments was performed for 45 cycles, followed by 25 cycles of the second round using a nested primer pair that targeted a specific segment in the diluted previous amplicons. The multiplex PCR reaction mixture (10 μl) consisted of 1× PCR buffer for Ex Taq HS (Takara), 0.25 mM dNTP (Takara), 0.6 μM of each primer (F15969, R16255, F16186, R16455, F16565, R289, F154, and R421), 2.5 mM MgCl2 (Roche), 0.5 U Ex Taq HS, 0.5× SYBR Green I (Sigma), 1.5 mg/ml BSA (Life Technologies), PCR-grade water (Roche), and 3 μl of aDNA extract. The multiplex PCR cycling conditions were as follows: 1 cycle of 2 min at 94°C; 45 cycles of 10 s at 94°C, 60 s at 57°C, and 80 s at 72°C with a single fluorescence acquisition at the end of the extension step; and one melting curve cycle of 10 s at 97°C, 120 s at 72°C, and a temperature increase from 65 to 97°C with a temperature transition rate of 0.1°C/s with continuous fluorescence acquisition at the rate of 5 readings/°C. The PCR products were diluted 250 times with PCR-grade water (Roche). For the second round of real-time PCR (nested PCR), 2 μl of the diluted products was used for the template, and nested primers for each target segment were added to the four separate reactions. The nested PCR reaction mixture (30 μl) consisted of 1× PCR buffer for Ex Taq HS (Takara), 0.2 mM dNTP (Takara), 0.2 μM of each nested primer (Table 2), 2 mM MgCl2 (Roche), 0.25 U Ex Taq HS, 0.5× SYBR Green I (Sigma), 0.5 mg/ml BSA (Life Technologies), PCR-grade water (Roche), and 2 μl of the diluted first-round product. The nested PCR cycling conditions were as follows: 1 cycle of 2 min at 95°C, 25 cycles of 10 s at 95°C, 10 s at 60°C, and 30 s at 72°C with a single fluorescence acquisition at the end of the extension step; and one melting curve cycle with the same conditions as described above. For the amplification of coding region segments, the method was the same as above except for the primers. Amplicon identification and purity were estimated based on melting curve analysis. The PCR products were purified. Sanger dideoxy sequencing was conducted with the nested primers bidirectionally by Bionics Incorporation, Korea. The sequences were read using SeqMan version 5.03 (DNASTAR) from the bidirectional sequence electropherograms of each DNA fragment. The sequences were uploaded to GenBank and assigned accession numbers.
Haplogroup predictionmtDNA haplogroups of aDNA samples were determined by using the Haplotracker application. First, the four control region sequences of each sample were used for the haplogroup prediction. The sequences were entered with the control buttons pressed and submitted to ‘Track by Sequence.’ Based on the most recent common ancestor (MRCA) and top-ranked haplogroups at the tracking results, superhaplogroups were predicted for subgroup level 0 or 1. To confirm them, their haplogroup-specific variants in the coding region were searched using the ‘Haplogroup database.’ Primers for the PCR amplification and sequencing of the segments located at the variant positions were designed as above. These coding region sequences and previous control region sequences were entered together into ‘Track by Sequence’ to confirm superhaplogroups of samples.
Quantitative and quality assessments of aDNA mtDNA moleculesReal-time PCR assay was performed using a LightCycler 96 (Roche) system to quantify the number of mtDNA molecules and assess their quality in the aDNA extracts as follows. A 287 bp segment of HV1 mtDNA region was assayed using the primer sequences shown in Table 2. The PCR reaction mixture was the same as the nested PCR reaction mixture described above, except for the reaction volume (10 μl), the primers, and 2 μl of aDNA extracts. Thermal-cycling conditions were the same as the above-described multiplex PCR cycling conditions, except for the extension time (60 s) and the number of cycles (40 cycles). Product specificity was controlled using the melting-curve analysis of the LightCycler 96 software (version 1.1.0.1320) (Roche). Ten-fold serial dilutions from 6 × 106 to 60 copies) of purified 278 bp PCR product quantified by using an Epoch spectrophotometer system (BioTek instruments) were used as templates to generate standard curves. The real-time PCR experiment was performed in duplicate with two ‘non-template controls.’ Analysis of the data was performed using LightCycler 96 software (version 1.1.0.1320) (Roche) to generate standard curves to calculate the DNA amount from each unknown sample. Standard curves showing the correlation coefficient of the trendline higher than 0.95 were used. The degree of PCR inhibition was evaluated based on the quantification cycle (Cq) and end-cycle fluorescence (ECFL) compared with the results of a human reference DNA (K562, Promega). The Cq was used for the quantification of the template DNA in quantitative PCR (qPCR). The Cq increased as the template quantity decreased or PCR inhibition occurred. The ECFL was defined as the fluorescence at the final cycle of amplification, and reflected the amount of the final PCR product. The ECFL decreased as PCR inhibition occurred (Kim et al., 2013).
We previously developed an improved ancient DNA extraction and purification using ion-exchange columns (Kim et al., 2008). This method has been effective for bone samples that were severely decayed and contaminated with dark-soil-origin substances. These samples often had such a low copy aDNA that large amounts of bone materials were required to recover sufficient amounts of aDNA templates for PCR. For these samples, silica-column-based methods, such as the QIAquick PCR Purification Kit, were not able to completely remove PCR inhibitors from the extracted DNA still showing dark in the purified elutes and consistent PCR failures. Only ion-exchange column methods were able to remove the dark color from the DNA extracts and the subsequent PCR worked successfully. However, this method required extended extraction time, and was laborious and hazardous due to the use of the extremely corrosive reagents. Instead, we developed an efficient, safe, and less laborious method based on previously reported method (Rohland and Hofreiter, 2007) but with slight modifications. The chief modification was that we used Qiaquick PCR purification at the final step of the aDNA elution from the silica particles. This step contributed to concentrating diluted aDNA elutes and removing residual silica particles, which are strong PCR inhibitors. Real-time quantitative PCR (Figure S1, Table S1) performed for mtDNA quantification and quality assessment of the aDNA extracts showed that the standard curves for quantification were acceptable with a correlation coefficient of 0.99. The amount of mtDNA in the aDNA extracts was estimated to average 12000 (1300–35000) copies per gram of bone powders. The amount of mtDNA in the aDNA templates used for PCR was estimated to be an average of 464 (52–1400) mtDNA copies per microliter of aDNA extract sufficient for PCR amplification (Figure S1, Table S1). The ECFL values of the aDNA extracts used to evaluate the degree of PCR inhibition were measured to be 3.1 ± 0.6 (1.8 ± 0.1–3.9 ± 0.0) on average. Compared with the ECFL values of control human DNA dilutions (4.1 ± 0.0 or 4.0 ± 0.1), most aDNA extracts were found to be slightly inhibited during PCR. One MNW1 extract exhibited the lowest ECFL (1.8 ± 0.1), indicating the highest degree of inhibition among the samples. However, since this sample has a relatively high number of mtDNA copies (1400 copies), using a dilute extract as a template may help to lower the degree of inhibition and increase the amount of amplicons after completion of PCR. It will also save DNA samples. We used Haplotracker to find the most conserved regions for primer position to secure sufficient primer hybridization, and thereby achieve successful PCR regardless of sample haplogroups by avoiding the area of primer-template mismatching sequences among all the target haplogroups. We designed multiplex PCR primers to simultaneously amplify the mtDNA control region segments, HV1 and HV2, to generate amplicons smaller than 300 bp. This multiplex PCR was efficient at obtaining multiple target sequences simultaneously in a small reaction volume with a reduced amount of aDNA templates. This was desirable for saving valuable aDNA samples. We also designed the primers to avoid sequences after poly-cytosine stretch that are commonly found in the HV1 and HV2 regions and interfere with obtaining clean sequences. We established optimal conditions for reaction mixtures and ran cycles for the multiplex real-time PCR with aDNA samples. We included Ex Taq HS, which has previously been demonstrated to be excellent for PCR amplification from the skeleton origin aDNA (Kim et al., 2015). The concentrations of reaction materials such as MgCl2, deoxynucleotide triphosphates (dNTP), polymerase, bovine serum albumin (BSA), and primers were optimally adjusted for multiplex PCR of aDNA samples. The Real-time PCR method allowed for a quick screening of successful production of target DNA segments without the need of electrophoresis, reducing labor and time.
We successfully obtained clean sequences from all target mtDNA segments of all aDNA samples tested in this study, HV1 and HV2 segments of control region, and coding region segments (Table S2). The sequences were obtained from two different bone fragments or teeth of each sample and bidirectionally from direct PCR sequencing. There were no inconsistent results between the four sequences of the target DNA segment in each of all aDNA samples in this study. No heteroplasmy sequences were found, except for sample MNE2. This sample had a length heteroplasmy (LH) between nucleotide position (np) 16184 and 16193, according to the revised Cambridge reference sequence (rCRS), which has been reported to be one of the commonly found LHs in mtDNA sequences of human individuals (Bendall and Sykes, 1995; Berger et al., 2011). One extract was sequenced with 10 cytosines, the other with 11 cytosines. However, this LH did not affect haplogroup classification. We treated the LH sequence as nine cytosines that appeared consistently clean among the duplicated samples. We predicted super-haplogroups (level 2) of the aDNA samples using Haplotracker with the control region sequence fragments HV1 and HV2 (Table 3). Haplotracker is an mtDNA haplogrouping web application developed based on the Phylotree build 17 (https://www.phylotree.org) (van Oven, 2015), the tree of global mtDNA variation. Haplotracker with control region segment sequences predicted a single haplogroup for five samples (MNE1, MNE2, MNE3, MNW1, and MNW4) in rank group no. 1. It predicted five haplogroups for MNW3 but haplogroup G1 at the top rank. None of the predicted haplogroups were identical. To confirm the haplogroups, haplogroup-determining variants in the coding region for each predicted haplogroup were sought using the Haplotracker Database tool. The coding region variants used for haplogroup determination are shown in Table 3. One or two coding-region variants were used for each haplogroup. All haplogroups confirmed with the haplogroup-determining variants were identical to those predicted using control region sequences. These results suggested that the potential utility of control region sequence fragments using Haplotracker in mtDNA haplogroup prediction.
| Samples | Control region | Coding region | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| HV1 | HV2 | Predicted haplogroupa | Range | Variantsb | Confirmed haplogroup | |||||
| Range | Variants | Range | Variants | |||||||
| MNE1 | 16003–16196 | 35– | 73 263 | G2 (1) | 13507–13605 | 13563 (G2) | G2 (1) | |||
| 16217–16413 | 16223 16227 16234 16278 | 365 | 309.1C | 5051–5167 | 5108 (G) | |||||
| 16362 | 315.1C | |||||||||
| MNE2 | 16039–16192 | 16189 | 35– | 73 152 | A+152 (1) | 1690–1815 | 1736 (A) | A+152 (1) | ||
| 16217–16413 | 16223 16262 16290 16319 | 392 | 235 263 | |||||||
| 315.1C | ||||||||||
| MNE3 | 16003–16196 | 47– | 73 200 | D4 (1) | 8408–8503 | 8414 (D4) | D4 (1) | |||
| 16217–16413 | 16223 16301 16319 16362 | 392 | 239 263 | |||||||
| 297 | ||||||||||
| 315.1C | ||||||||||
| MNW1 | 16003–16196 | 35– | 73 189 | W+194 (1) | 15728–15901 | 15784 15884C | W+194 (1) | |||
| 16217–16413 | 16223 16292 | 365 | 194 195 | (W) | ||||||
| 204 207 | ||||||||||
| 263 | ||||||||||
| 309.1C | ||||||||||
| 315.1C | ||||||||||
| MNW3 | 16003–16196 | 35– | 73 150 | G1 (1) | 15280–15390 | 15301 | G1 (1) | |||
| 16217–16413 | 16223 16325 | 365 | 263 | L3e’i’k’x L3e L3i | 15323 (G1) | |||||
| 315.1C | M62'68 (5) | 15326 | ||||||||
| 5051–5167 | 5108 (G) | |||||||||
| MNW4 | 16003–16196 | 16129 | 35– | 73 195 | C4 (1) | 15133–15207 | 15204 (C4) | C4 (1) | ||
| 16217–16413 | 16223 16298 16327 | 365 | 249d | 3493–3598 | 3552A (C) | |||||
| 263 | ||||||||||
| 309.1C | ||||||||||
| 315.1C | ||||||||||
HV1, hypervariable region 1; HV2, hypervariable region 2.
Variants were written according to the notation in Phylotree Build 17 (https://www.phylotree.org/). Haplogroup-characteristic variants are shown in bold.
a Haplogroups were predicted at super-haplogroup level 2 using Haplotracker (https://haplotracker.cau.ac.kr); ranks are in parentheses.
b Variant’s specific haplogroup is shown in parentheses.
Xiongnu empire cemeteries of the elites were located throughout Mongolia. The Xiongnu cemetery in Duurlig Nars was first discovered by a group directed by D. Tseveeendorj in 1974. It is considered the biggest Xiongnu cemetery in eastern Mongolia (National Museum of Korea et al., 2011). Tombs in Duurlig Nars have been excavated with various gold ornaments; therefore, it has been considered as the cemetery of high-class society showing the social stratification of the Xiongnu empire (Kim et al., 2010; National Museum of Korea et al., 2014).
With the excavation of three satellite tombs from the 2000-year-old Xiongnu cemetery in Duurlig Nars in 2006–2007, we investigated the kinship between the three ancient human skeletons from the tombs by molecular genetics study. We analyzed maternal mtDNA, paternal Y-chromosome single nucleotide polymorphisms (Y-SNP), and biparental autosomal short tandem repeats (STR). The two ancient humans were one female with the mtDNA haplogroup D4 and one male with the Y-SNP haplogroup C3 and mtDNA haplogroup D4. One ancient human was a male with the maternal U2e1 and paternal R1a1 haplogroup. A male of Indo-European lineage (R1a1) was present in the Xiongnu region of Mongolia 2000 years ago. Moreover, there was no immediate kinship among the three aDNA samples, with different mtDNA and Y-SNP haplogroup and sequences (Kim et al., 2010).
Main tombs with seven satellite tombs were excavated in 2007 in the Tahiltin-hotgor cemetery known to be Xiongnu elite tomb complexes and located in Manhan sum of Hovd Aimag, Mongolia. A molecular genetic study was not performed at this archaeologic site (Jones and Joseph, 2008; Miller et al., 2008). At Arkhangai Aimag, Mongolia, the Gol Mod 2 Xiongnu Elite cemetery was discovered in Ondor-Ulaan sum. Around the main standard Xiongnu tomb, 30 satellite burials were found; however, no molecular genetic analysis of the satellite tombs was reported (Erdenebaatar et al., 2011). One main tomb with three satellite tombs from the Xiongnu period was excavated in 2017 by a Korea–Mongolia joint team. However, there was no further genetic investigation of human skeletons from the four tombs (Oh et al., 2018). A study was performed in the 2000-year-old Egyin Gol Xiongnu cemetery in northern Mongolia using biparental, paternal, and maternal genetic systems. Kinship was constructed using the molecular genetic analysis of human skeletons excavated from the cemetery (Keyser-Tracqui et al., 2003). However, the tombs in the cemetery seemed not to show satellite sites. Instead, they showed aggregated sites of family or relatives of ordinary peoples. The Tamir Ulaan Khoshuu cemetery in the Arkhangai Aimag in central Mongolia, was excavated by a French–Mongolian joint research team between 2013 and 2018. Combined analysis of nuclear and whole mitochondrial DNA was performed on the skeletal remains of 52 individuals excavated from 48 stone ring burials. The mitochondrial DNA data showed 11 different main haplogroups of European (H, J, T, U, X) and Asian (A, B, C, D, F, G) origin. The at least five Y haplogroups were revealed (R, Q, N, J and G). A kinship of a large family was constructed spanning five generations with two male and six female haplotypes (Keyser et al., 2021).
Xiongnu territory is overlapped with Mongolia; therefore, many Xiongnu cemeteries are present in Mongolia. Terrace tombs are generally large and considered as elite with many valuable funeral goods found within them despite pillaging. Rectangular terrace tombs with satellite graves—small circular tombs found around the main large terrace tombs—have been frequently found throughout Xiongnu cemeteries. The individuals in the satellite tombs can be members of a family, including concubines, relatives, loyal retainers, interior subordinates, bodyguards, servants, slaves, or sacrifices. To date, the relationship between the bodies in satellite tombs and those in the main terrace tomb have not been elucidated.
To the best of our knowledge, there has been no molecular genetic analysis of ancient human subjects using mtDNA haplogroups in satellite tombs in Xiongnu elite cemeteries. This is the first study to report the relationship between satellite tombs around the main Xiongnu elite tomb. The mitochondrial DNA haplotypes showed six different haplogroups (A, C4, D4, G1, G2, and W) in the human skeletons of six satellite tombs. The maternal origins are all different. All the mtDNA haplotypes were Asian except the W. W haplogroup in ancient human specimens has been found in Germany, Spain, and Kazakhstan (Olivieri et al., 2013). The first expansion of haplogroup W occurred in South Asia (Kivisild et al., 2003) and the secondary expansion spanned from Near East Asia to western India (Endicott et al., 2007).
Double graves in a 2000-year-old Xiongnu necropolis could be a sacrificial burial (Keyser-Tracqui et al., 2003). Sacrificial burial was also suggested for an ancient tomb containing a woman, where the hyoid bone found in an offering box (Murail et al., 2000). In these satellite tombs, five ancient human skeletons were studied anthropologically (Kim, 2014). Three skeletons (T1, E1, W1) were estimated as males in their sixties, and two skeletons as females in their twenties (E2, E3). The E1 and E2 tombs were located side by side; therefore, the E2 tomb could be a sacrificial burial.
In summary, we have developed an efficient and safe DNA extraction method which is less laborious than the previously reported method using ion-exchange columns. The method is based on the report by Rohland and Hofreiter (2007) but with some modifications. We designed multiplex PCR to simultaneously amplify the mtDNA control region segments, HV1 and HV2. This was efficient at obtaining multiple target sequences and desirable for saving valuable aDNA samples. This is the first study to report the relationship between ancient human skeletons using molecular genetics in the satellite tombs throughout the Xiongnu empire. The various maternal lineages in the satellite tombs around the Xiongnu elite tomb suggested that Xiongnu society was permissive towards others outside its bloodline.
Mitochondrial analysis was performed on six Xiongnu satellite tombs at the 2000-year-old Xiongnu elite cemetery in Duurlig Nars, Mongolia. This is the first study to report the genetic relationship between individuals in satellite tombs. The different mitochondrial haplotypes of ancient human skeletons from six satellite tombs suggested that there was no maternal kin relationship between the six individuals. One sacrificial tomb was proposed, which was a female in her twenties. Further studies should use Y haplotypes, Y chromosomal STRs, autosomal STRs, and phenotype and biogeography ancestry–informative DNA.
The data used in this study are provided in the figures, tables, and supplementary materials.
All new DNA sequences found in this study have been submitted to GenBank. GenBank accession numbers of the sequences are presented in Table S2. The deposition is under processing.
The authors have no conflicts of interest to declare.
This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF2018R1D1A1B0505025214).
We thank the officials of National Museum of Korea for providing us with ancient Mongolian human samples.
