In total, bloodstains within FTA™ cards were collected from 91 unrelated healthy individuals of the Kirgiz group in the Northwest region of China. The detailed sampling locations of all 49 populations including Kirgiz group and 48 reference populations were shown on the world map in Figure 1. All participants were informed in details of the content and purpose of the present study, and all signed the written informed consents. This study was approved by the Ethics Committee of Xi’an Jiaotong University Health Science Center (Approval Number: XJTULAC201) and complied with the ethical principle of the World Medical Association Declaration of Helsinki (World Medical Association, 2013). According to the questionnaire survey, the selected donors of random sampling were all local residents and their immediate family members (parents, grandparents, and maternal grandparents) were all Kirgiz individuals. To protect the privacy of volunteers, all samples were numbered and anonymized during the whole experiments.

Genome mtDNA of the 91 bloodstain samples was extracted with the PrepFiler BTA™ forensic DNA extraction kit (Thermo Fisher Scientific, MA, United States). The magnetic bead method was used for DNA extraction by the above kit. The extraction method selected in this study combined uniquely structured magnetic particles and the optimized multi-component chemistry surface, to provide efficient DNA binding capacity and high DNA recovery rate. The extracted genome DNA was quantified with the Qubit™ dsDNA HS assay kit (Thermo Fisher Scientific) at the Qubit™ four Fluorometer (Thermo Fisher Scientific). The mtDNA control region sequences of 91 unrelated individuals and positive control samples (007 and 9947A) were amplified using the Precision ID mtDNA Control Region Panel (Thermo Fisher Scientific). The libraries were prepared with Precision ID DL8 Kit (Thermo Fisher Scientific) on Ion Chef™ System (Thermo Fisher Scientific). The mtDNA control region was amplificated with 14 primer pairs in two primer pools. The libraries were quantified using the Ion Library TaqMan™ Quantitation Kit (Thermo Fisher Scientific) on the 7,500 Real-Time PCR System (Thermo Fisher Scientific). The sequencing templates were prepared using with the Ion S5™ Precision ID Chef & Sequencing Kit (Thermo Fisher Scientific), and were loaded to the Ion 530™ Chip (Thermo Fisher Scientific) on the Ion Chef™ System (Thermo Fisher Scientific). The sequencing runs were accomplished using the Ion S5™ Precision ID Chef & Sequencing Kit (Thermo Fisher Scientific) on Ion S5™ XL System (Thermo Fisher Scientific). All the above experimental procedures were performed strictly according to the manufacturer’s recommendations.

The raw sequencing data were analyzed using Converge ™ Software V2.2 (Thermo Fisher Scientific) of MPS mtDNA module with HID Genotyper 2.2 of default parameters. After aligning to the Revised Cambridge Reference Sequence (rCRS) (Andrews et al., 1999), all the variations were called for further scoring. After the mutations were verified, they were mapped to nodes in the Phylotree 17 (http://www.phylotree.org/tree/index.htm) of human mitochondrial DNA lineages, further allocated the haplotypes to the closest haplogroups, which were submitted to the EMPOP database for quality control (Parson and Dür, 2007). The haplotypes were annotated according to the recommendations of International Society of Forensic Genetics with the nomenclature rules of mtDNA typing (Parson et al., 2014). The analytical threshold at all 1,122 sequences was 20×, and the thresholds of confirmed mutation, point heteroplasmy, insertion and deletion (length heteroplasmy) were 96%, 10%, 20%, and 30% of the total reads with each variant, respectively.

A total of 5,795 individuals from 48 reference populations distributed all over the world were chosen for further investigation of the genetic relationships among the studied Kirgiz ethnic group and the reference populations based on the mtDNA control region sequences. The sample size of each selected population was more than 30 individuals. The 48 reference populations included seven populations from Africa, five populations from North America, six populations from South America, nine populations from Europe, four populations from Central Asia, eight populations from East Asia, four populations from South Asia, two populations from Southeast Asia and four populations from West Asia. The detailed information of the Kirgiz group and these 48 reference populations was showed in Supplementary Table S1.

To evaluate the sequencing performance, we calculated several parameters including the true allelic ratio, noise ratio, ratio of two sequencing directions and the mean depth of each sample, the average depth of each amplicon, and the depth ratio between two primer pools of each sample. The true allele ratios were calculated by dividing the read depth of the true allele into total read depth of each sample, and the remaining undefined reads were defined as noise. The sequencing ratios were calculated by dividing the coverage of 5′ reads into the 3’ reads. The mean depth per sample and average depth of each amplicon were also calculated directly. The depth ratios were by dividing the lowest read depth into the greatest read depth of each sample.

The haplogroup frequencies based on the mtDNA control region sequences of 91 Kirgizs were calculated directly. The phylogenetic tree for 91 Kirgiz mtDNA control region sequences was performed using the online tool of HaploGrep 2 v2.0 (https://haplogrep.i-med.ac.at/app/index.html) under the Kulczynski measure. The example case of U1a2 with detailed calculation formula could be found at the website https://haplogrep.i-med.ac.at/2018/06/21/explaining-the-formula/. These mutations including 309.1C(C), 315.1C, 515-522 InDel of AC, A16182C, A16183C, 16193.1C(C), C16519T and T16519C were excluded when the phylogenetic tree was reconstructed, and the tree was showed in Supplementary Figure S1.

DNA Sequence Polymorphism v6 (DnaSP v6) (Librado and Rozas, 2009) was used to calculate the statistical parameters of the 91 mtDNA control region sequences, which included the number of polymorphic site (S), haplotype diversity (Hd), nucleotide diversity (Pi), and the average number of nucleotide difference (k). Analysis of molecular variation (AMOVA), neutrality tests (Tajima’s D and Fu’s Fs tests), and mismatch distribution (the sum of square deviation, SSD, and Harpending’s raggedness index, R) were performed with mtDNA control region sequences using Arlequin 3.5.1.3, simultaneously (Excoffier and Lischer, 2010). The Bayesian phylogenetic inference and Bayesian Skyline Plot (BSP) of Kirgiz group were performed using BEAST 2.6.5 software to deduce the time when the Kirgiz group expansion occurred (Bouckaert et al., 2019). The nucleotide substitution model of TN93 and gamma site model with the substitution rate of 1.57E-8 of the single nucleotide substitution per year were set (Soares et al., 2009), and the strict clock model was selected. The pairwise divergence time for Homo sapiens and Homo neanderthalensis was chosen to calibrate the divergence time, as a median time of 0.55 (±0.054) million years ago (MYA) (Soares et al., 2009). The adaptive Monte Carlo Markov chain (MCMC) approach was used to evaluate the molecular evolution of the mtDNA control region sequences, the MCMC chain length was 1,000,000, and auto optimize displayed for 10,000 steps. Bayesian Skyline Plot (BSP) reconstruction was performed using Tracer 1.7.0 (Rambaut et al., 2018) and the phylogenetic tree was annotated by FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/). The detailed geographical coordinates and the macrohaplogroup compositions of all the populations were visualized using Tableau Desktop v 2021.1. Heatmap was constructed using the ‘heatmap’ package of R statistical software based on the pairwise F _ST values between the Kirgiz group and 48 reference populations (https://www.r-project.org/). A neighbor-joining (NJ) tree was built using Molecular Evolutionary Genetics Analysis (Mega) v10.2.1 based on the pairwise F _ST values among 49 populations. And the NJ tree was displayed and annotated using the online tool iTOL v 6.1.2 (https://itol.embl.de/itol.cgi). The two-dimensional (2D) and three-dimensional (3D) principle component analyses (PCA) were constructed using the dimensionality reduction analysis module of SPSS Statistics v.19.0 based on the haplogroup frequencies generated from mtDNA control regions of the Kirgiz group and 48 reference populations. The haplogroup network calculations were performed using Network v.10.2.0.0 of the median joining module based on the mtDNA control region sequences from the Kirgiz and the 48 reference populations, finally, of which 1,403 individuals from the 48 reference populations sharing the same or similar haplogroups with the target Kirgiz individuals were selected to participate in the network relationship evaluations. The following network plots were visualized with the network publisher (https://www.fluxus-engineering.com/nwpub.htm).

All the mtDNA control region sequences from the 91 Kirgiz individuals and the positive control samples (9947A and 007) were successfully generated from four sequencing runs on the Ion S5™ XL system. The true allelic ratios and noise ratios were shown in Figures 2A,B. Even when outliers were taken into consideration, the true alleles of all tested samples covered more than 96% of sequencing reads, while the noise was less than 4%. The balance ratios of all the amplicons from two directions in the 91 samples were shown in Figure 2C. The balance parameters of most samples were around 0.5, which indicated the relatively good performance of the bi-directional sequencing. In Figure 2D, the blue and green bars represented the average depths of two primer pools with 14 amplicons when the control region of mtDNA was amplified using the imbrication strategy. The mean depth per individual (the first red boxplot in Figure 2D) was 12989 ± 10865 × (mean ± standard deviation). The mean depths of the first pool ranged from 4,598 × (P1A6) to 24535 × (P1A3) which were shown in the blue boxplot in Figure 2D. The mean depths of the second pool ranged from 3,632 × (P2A5) to 33856 × (P2A1), as shown in the green boxplot in Figure 2D. The average depth ratio of two pools (the purple boxplot in Figure 2D) was 68.65%.

A total of 66 haplogroups and 81 haplotypes were identified from 91 Kirgiz unrelated individuals based on mtDNA control region sequences. The variants (including single point mutation, InDel, and heteroplasmy were marked in green, grey and yellow, respectively), haplotypes and allocated haplogroups of 91 Kirgiz individuals were shown in Supplementary Table S3. The mtDNA haplogroup distributions of the Kirgiz group were shown in Figure 3A. The haplogroups of 91 Kirgizs comprised of 69.23% East Eurasian haplogroups and 30.77% Western Eurasian haplogroups. The macrohaplogroup M (2.20% of Z1, 4.40% of M8, 5.49% of M, 4.40% of G, 2.20% of D5, 10.99 of D4, 1.10% of D1, 2.20% of C7, 3.30% of C5, and 9.89% of C4), macrohaplogroup R (5.49% of F1, 1.10% of B5, 5.49% of B4, and 3.30% R) and macrohaplogroup N (1.10% of N9, 6.59% of A) commonly observed in East Eurasia populations were detected in the Kirgiz group, respectively. Additionally, the Western Eurasian haplogroups including H (7.69%), HV (1.10%), J1 (2.20%), N1 (2.20%), N1 (2.20%), T (6.59%), U (4.40%), W6 (3.30%), X2 (2.20%), and R0 (1.10%) were also observed in the Kirgiz group.

DNA polymorphisms of the studied Kirgiz group were evaluated, and the results of indexes were presented as follow, number of polymorphic site (S: 107), number of haplotype (h: 81), haplotype diversity (Hd: 0.997), nucleotide diversity (Pi: 0.00842) and the average number of nucleotide difference (k: 9.017), these parameters indicated the relatively high genetic diversities of mtDNA control region sequences was in Kirgiz group. The Tajima’s D test (-1.950, p-value < 0.05) and Fu’s FS test (-25.247, p-value < 0.05) of mtDNA variants for the 91 Kirgiz individuals were both calculated to be relatively large negative values, and the results might indicate that the variants were significantly deviated from neutral mutations.

The mismatch distribution graph shown in Supplementary Figure S2 revealed a single peak, which indicated the Kirgiz group might occur the population expansion event according to the previous research method (Slatkin and Hudson, 1991). The SSD and Raggedness indexes of two models in the Kirgiz group were 0.00182 (SSD, p-value = 0.60000), 0.00689 (Raggedness index, p-value = 0.73000) for the sudden expansion model, 0.00216 (SSD, p-value = 0.57000) and 0.00689 (Raggedness index, p-value = 0.86000) for the spatial expansion model, respectively. The median network calculation was performed based on the mtDNA control region sequences of the 91 Kirgiz individuals, and the network plot of 91 Kirgiz individuals was shown in Figure 3B. Small branches belonged to the same haplogroup clustered together and separated from the other different haplogroup branches. Furthermore, as shown in Figure 4, Bayesian Skyline Plot (BSP) analysis was conducted based on the ancestral theory to quantify the evolutionary background of population size and history. The expansion time with the largest slope of the BSP abscissa was directly read. And the BSP abscissa showed that the Kirgiz group had grown at about 53.41 kya (Figure 4A, the abscissa corresponding to the dash line). The annotated phylogenetic tree (Figure 4B) showed the estimated divergence time between Homo sapiens and Homo neanderthalensis was 0.5413 MYG, which was very close to the expected time of 0.55 MYG.

As for an important genetic marker for population evolution analysis and ancestral information inference, mtDNA control region sequences originated from 91 Kirgiz individuals were combined with other 5,795 mtDNA control region sequences from those individuals of 48 reference populations, which were used to analyze genetic relationships between Kirgiz group and these reference populations. The detailed geographical coordinates and the macrohaplogroup compositions of the Kirgiz group and 48 reference populations were shown in Figure 1 and Supplementary Table S4. As shown in Figure 1, with the advancement towards Africa, Europe, Asia, and America, the proportion of African-specific macrohaplogroup L gradually decreased and disappeared, and the proportions of macrohaplogroups M, N and R gradually increased. In macrohaplogroups N and R in this figure, the Western Eurasian specific haplogroups had been excluded. The 92.77% haplogroups of Caucasian population from North America (CAU) belonged to the Western Eurasian haplogroup because of the European origin.

Pairwise F _ST values were calculated based on the haplogroup frequencies of mtDNA control region sequences to obtain insight into the genetic affinities of the Kirgiz group and 48 reference populations. A population clustering heatmap of the pairwise F _ST values was plotted and shown in Figure 5A. In addition, the NJ tree of these 49 populations was displayed in Figure 5B. In the NJ tree, the Kirgiz group gathered together with other East Asia populations and formed a larger cluster. As a whole, the NJ tree indicated that most geographically adjacent populations gathered together. However, some populations were outliers which might be due to the population histories or the deviations caused by the different methodologies. The pairwise F _ST values were shown in Supplementary Table S6. The pairwise F _ST values between the Kirgiz group and other reference populations ranged from 0.01142 to 0.1499, and the median was 0.02836. The largest three F _ST values with the Kirgiz group were identified in three populations including Peruvian in Lima (PEL), Valles (VAL) and Colombian in Medellin (CLM). Among Chinese nine populations, the smallest three F _ST values were observed between Kirgiz and Beijing Han Chinese (CHB), East Pamiri Kyrgyz (EPK) and Southern Han Chinese (CHS).

PCA was conducted to evaluate the genetic affinities among Kirgiz group and reference populations based on the haplogroup frequencies in terms of maternal inheritance, and the haplogroup distributions of Kirgiz group and 48 reference populations were shown in Supplementary Table S5. And the scree PCA plot was shown in Supplementary Figure S3, which represented the principal components and percentages of explained variances. The first three principal components could explain 38.82% of the total variation among these 49 populations, while the first principal component explained 17.87% variation, the second explained 10.66%, and the third explained 10.29%. Based on the first three principal components, a 3D-PCA (PC1, PC2 and PC3) and a 2D-PCA (PC1 and PC2) were constructed and presented in Figure 6. The populations from nine different continental regions were marked with nine different colors. The PCA plots in Supplementary Figure S4 showed the 2D-PCA plots (PC1 and PC3, PC2 and PC3). As displayed in the PCA plots, the East Asia populations could be separated in PC1, the American populations could be separated in PC2, and the African populations could be separated in PC3.

To evaluate the network relationships between Kirgiz and other reference populations, all the 91 Kirgiz mtDNA control region sequences were compared with the individual sequences from 48 reference populations, and six median network graphs were finally generated. In total, 460 different individual sequences of macrohaplogroup M were distributed into seven different haplogroups, as 62 of D4, 36 of M1, 216 of M, 44 of G, 63 of M8, eight of D1, and 31 of D5, respectively. The median network diagrams of macrohaplogroup M with the largest proportion of Kirgiz group were labeled by different colors for different geographical locations, and different haplogroups were shown in Figure 7 and Supplementary Figure S5, respectively. In Figure 7, the red arrows indicated the 42 Kirgiz individual sequences were allocated into macrohaplogroup M. There were 58 individuals from eight populations from East Asia, one population from South Asia, two populations from Southeast Asia, one population from Europe, three populations from South America and two populations from North America which had homogenous haplogroup M with three Kirgiz individuals. Two East Asia populations were observed the direct link with Kirgiz group. Haplotype sequences apart from one median vector were observed between the Kirgiz individuals and one Southeast Asian, three East Asian individuals.

The commonly observed Western Eurasian haplogroup including the H, T, U, W, X2, J1, N1, HV and R0 was the second-largest haplogroup in the Kirgiz group. These 27 Kirgiz individuals were allocated into this category. The Western Eurasian haplogroup of 27 Kirgiz individuals was allocated into nine above-mentioned haplogroups and nine different network plots labeled with colors for different continents (Figure 8) and different haplogroups (Supplementary Figure S6), respectively. There were 582 mtDNA control region sequences including in these nine median network relationships: 80 of H, 77 of T, 65 of U, 64 of W, 42 of X2, 97 of J1, 74 of N1, 47 of HV, and 36 of R0. Compared with the Kirgiz group, the directly linked individuals were observed in five European populations including Toscani in Italy (TSI), Northern Greece (NOG), Finnish in Finland (FIN), Iberian populations in Spain (IBS) and Moldavia (MOL), two West Asia populations i.e. Kuwaiti (KUW) and Jordan (JOR), one South Asia group Bengali in Bangladesh (BEB), one South America population Lebanon (LEB) and one African population Arab (ARA). Individuals apart from one median vector with Kirgizs were mostly observed in European populations, followed by the Central Asia and West Asia populations.

The 14 individuals of Kirgiz group were allocated into the macrohaplogroup R. Combining with 263 individuals from reference populations, 277 mtDNA control region sequences were distributed into four different haplogroups. And these four different network plots were shown in Supplementary Figure S7. These four different networks were marked by different colors for different continents in Supplementary Figure S7A, whereas, the different haplogroups were labeled by different colors in Supplementary Figure S7B, respectively. The median networks of four haplogroups B4, B5, R, and F1 were reconstructed based on the 91, 82, 70, and 20 mtDNA control region sequences among 48 populations, respectively. One individual of Utah resident with Northern and Western European ancestry (CEU) shared the homogenous haplogroup with a Kirgiz individual. Some individuals from eight populations including seven populations of Asia and one of Europe linked directly with 14 Kirgiz individuals, respectively. Seven individuals from three South Asian populations and one Southeast Asian, one North American, one South American, and one African population were found to be apart from one median vector with 14 Kirgiz individuals, respectively.

The smallest macrohaplogroup proportion of 91 unrelated Kirgiz individuals was N macrohaplogroup. The 25 sequences from the reference populations and seven sequences from the Kirgiz group were used to draw the median network diagram of the N macrohaplogroup. And the diagram of each haplogroup was labeled by different continents and haplogroups with different colors (Supplementary Figure S8). Within the macrohaplogroup N, seven Kirgiz individuals clustered together and separated from the individuals of other reference populations.