A total of 938 (360♂♂, 578♀♀) host-seeking adult D. marginatus individuals were collected from 151 pinpoint locations across 31 sampling sites. This included 553 individuals from 102 locations within 19 sampling sites in CN and 256 individuals from 49 locations within 12 sampling sites in NE. All individuals were classified into the D. marginatus complex based on morphological criteria and identified as the D. marginatus morphotype. For the population genetics analysis, a total of 361 D. marginatus individuals (168♂♂, 193♀♀) were selected from 31 sampling sites, with each site contributing between 8 and 21 individuals. Additionally, to ensure broad ecological and geographical coverage, at least one tick was collected from each of the 150 locations. However, the primary focus of this study is not on individual locations, but rather on the genetic differentiation between sampling sites and the two major regional populations: CN (n = 226) and NE (n = 135). Given this approach, the sample size is considered sufficient for assessing regional population structure. Detailed information about the collected tick specimens included in the genetic analysis is provided in Table 1 and Table S2, and their geographical distribution is illustrated in Fig. 1.

All D. marginatus samples included in the population genetics were subjected to PCR amplification of the partial cox1 and the entire ITS2 gene, as outlined in the methodology. Positive amplicons of the expected sizes were successfully obtained for all samples. The sequences from these amplicons were assembled bidirectionally, and after eliminating low-quality regions at the beginning and end, the final cox1 sequence had an average length of 840 bp. For the ITS2 gene, a complete sequence of 1099 bp was obtained. Additionally, to be used as an outgroup and for phylogenetic analysis, the entire ITS2 gene region of D. raskemensis was sequenced. Although cox1 data for D. raskemensis (GenBank accession MT308586) were already available, ITS2 data were absent in the GenBank database. Thus, ITS2 sequencing was performed on the same specimen from our DNA bank, and the resulting sequence was deposited in GenBank under the accession number PP618825.

Population structure and demographic history based on the cox1 gene

The cox1 sequences of 361 characterized individuals were aligned as amino acids using the MUSCLE algorithm and trimmed to match the shortest sequence, resulting in a dataset of sequences with a total length of 824 bp. As the cox1 gene is protein-coding, the nucleotide sequences were also translated into amino acid sequences and analyzed for the presence of stop codons and numts. No stop codons or numts were detected in any of the sequences. Analyses were performed separately for each sampling site, as well as for the two regional groups and the overall population. However, it was determined that the results were more informative when analyzed on a regional basis, so most inferences were drawn from the analyses of the two regions and the entire population combined.

Central Anatolia (CN)

In the CN, 71 haplotypes were identified, with 67 of which (CX-CN1-67) were unique to this region, while four haplotypes (CX-CNNE1-4) were shared between regions. Among these, 49 individuals were represented by a single haplotype, while the remaining samples were grouped into haplotypes with multiple representatives. Neutrality tests revealed negative and statistically significant values for Tajima’s D and Fu and Li’s D and F (P < 0.02), while the Fu’s Fs value was also negative (Table 2). The mismatch distribution analysis displayed a unimodal pattern (Fig. S1).

Northeast Anatolia (NE)

In the NE, 64 haplotypes were identified, including 60 unique haplotypes (CX-NE1-60) that are specific to this region and four haplotypes (CX-CNNE1-4) shared between regions. Among these, 49 individuals were represented by a single haplotype, while the remaining individuals were grouped into haplotypes with multiple representatives. Neutrality tests showed that Tajima’s D, Fu and Li’s D and F values were negative and statistically significant (P < 0.02), and Fu’s Fs value was also negative (Table 2). The mismatch distribution analysis revealed a unimodal pattern (Fig. S1).

All regions (ALL)

When all study areas were considered as a single population, a total of 131 haplotypes were identified. Of these, 67 were specific to CN, 60 to NE, and four were shared between the two regions, accounting for all haplotypes detected in the study. Among the samples, 98 individuals were represented by a single haplotype, while the remaining samples were grouped into haplotypes with multiple individuals. Neutrality tests showed that Tajima’s D, Fu and Li’s D and F values were negative and statistically significant (P < 0.02), and Fu’s Fs value was also negative (Table 2). The mismatch distribution analysis exhibited a unimodal pattern (Fig. S1). Additionally, the results of the analysis of polymorphic regions are presented in Table S3, based on both the total population (all combined samples) and the two regional populations (CN and NE).

Genetic differentiation and population structure

To determine genetic variance and pairwise genetic differentiation (F_ST) at the cox1 gene level, populations were categorized into 31 sampling sites and two regional populations. AMOVA was performed locus by locus. In the dataset with 31 populations, genetic variation among populations accounted for 15.94%, while within-population variation was 84.07%, yielding an F_ST value of 0.15936 (P < 0.001). In the two regional populations, genetic variation among populations was 19.35%, while variation within populations was 80.64%, with an F_ST value of 0.19353 (P < 0.001) (Table S4). The most differentiated population, with the highest F_ST value (0.41205, P < 0.001), was found between sampling site L6 in CN and L16 in NE. A distance matrix and color plot generated from the F_ST values for the 31 populations are provided in the supplementary files (Fig. S2 and Table S5).

STRUCTURE analyses were conducted with simulations ranging from K = 2 to 10 groups (27 simulations in total, three simulations for each group). The most appropriate number of groups was determined to be K = 3 using pophelper (Fig. S3), and the population was evaluated based on three distinct ancestral groups. This analysis revealed significant population structure between the CN and NE populations of D. marginatus. The STRUCTURE graph indicated that red alleles were predominant in the CN population, while green alleles dominated in the NE population, with black alleles being rare and recessive in both regions (Fig. 2A). The correlation between genetic variation and geographic distance was assessed using the Mantel test on the dataset of 361 D. marginatus individuals, based on cox1 gene sequences and individual geographical coordinates. The test yielded a statistically significant result (r = 0.2785, P < 0.05). Additionally, the dataset, with geographical coordinates based on 31 sampling sites, was subjected to SAMOVA analysis. Simulations from K = 2 to 10 groups were performed, with K = 4 being identified as the most appropriate based on the F_CT value. In the K = 2 simulation, the entire population was divided into two separate groups: CN and NE populations. As the number of groups increased up to K = 7, separations were observed within the NE populations, while the CN population remained as a single group. At K = 7, LGr from CN was the first to separate. At K = 4, the CN population formed a single group, while the NE population was divided into three groups: Group 1: L17; Group 2: L18, L19, L20, L21, L22, L23, L24, L25, L26, L27; Group 3: L16 (Fig. 2B).

Haplotype distribution, Haplotype network, and Pairwise-Distance

The 131 characterized D. marginatus individuals were represented by 131 distinct haplotypes based on the cox1 gene. Among these, 67 haplotypes were unique to CN, 60 haplotypes were specific to NE, and four haplotypes were shared between the two regions. In CE, the most common haplotype, CX-CN3, was represented by 77 individuals, followed by CX-CN2 with 35 individuals, and other haplotypes with progressively fewer individuals. The CX-CN3 haplotype was present in all sampling sites in CN, while CX-CN2 was found in all but three sampling sites (L7, LA, and LGr).

In NE, the most common haplotype was CX-NE15, represented by 10 individuals, followed by CX-NE1 with six individuals. These haplotypes were distributed across specific sampling sites, such as CX-NE15, which was found in L18, L20, L21, L22, and L27, and CX-NE1, found in L16, L18, L21, and L22. The shared haplotypes, CX-CNNE1, CX-CNNE3, CX-CNNE2, and CX-CNNE4, exhibited varying distribution patterns, with CX-CNNE1 being the most abundant, found in eight sampling sites in CN and across all sampling sites in NE. Detailed information on the haplotypes and their distribution across the different sampling sites is provided in Table S6 and illustrated in Fig. 3.

A pairwise-distance analysis was performed using the 131 identified haplotypes. The most genetically distant haplotype was CX-NE53, represented by a single individual (Dm-2411) from sampling site L24 in NE. Intraspecific genetic variation ranged from 0.12% to 1.94%, with noticeable pairwise differentiation between region-specific haplotypes. Genetic differentiation was more pronounced within haplotypes from the NE (Table S7 and Fig. 3).

The 131 characterized haplotypes were analyzed for genetic similarity and haplotype specificity using BLAST and identification analyses in both the GenBank and BOLD databases. In the GenBank database, the haplotypes could be compared with sequences of equal size, while in the BOLD database, comparisons were made with sequences that were, on average, at least 200 bp shorter than our haplotypes. Therefore, similarity and uniqueness analyses were based on BLAST results. According to the BLAST analysis, none of the haplotypes were identical to any existing records (whether through full or close length comparisons), confirming that all 131 haplotypes were unique. The most similar sequences, with homologies ranging from 98.30% to 99.88% (comparison rate above 98%), were found to be from D. marginatus in Kazakhstan (MN907848 and OQ415364), Slovakia (MK905212), and China (NC_062069 and OM368304). While comparison rates in the BOLD database were generally low, identification analysis showed 99.38% to 100% similarity with D. marginatus records under the AAL1447 BOLD accession, primarily from Spain, Romania, Croatia, and Georgia (Table S8).

Population structure and demographic history based on the ITS2 gene

The complete ITS2 gene sequence (total length: 1099 bp) of 361 characterized D. marginatus individuals was aligned using the MAFFT algorithm. Comparisons within the dataset revealed no insertions or deletions among the sequences. Analyses were conducted separately for each sampling site, as well as at the levels of two regions and the entire population.

Central Anatolia (CN)

In the CN, 65 genotypes were identified, comprising 37 region-specific genotypes (IT-CN1-37) and 28 common genotypes (IT-CNNE1-28). Among these, 33 individuals were represented by a single genotype, while the remaining samples were grouped into genotypes shared by multiple individuals. Neutrality tests revealed mixed results: Tajima’s D was positive and statistically insignificant, Fu and Li’s D was positive but statistically significant (P < 0.05), Fu and Li’s F was negative and statistically insignificant, and Fu’s Fs was negative (Table S9). The mismatch distribution analysis exhibited a unimodal-like pattern (Fig. S4).

Northeast Anatolia (NE)

In the NE, 67 genotypes were identified, including 39 region-specific genotypes (IT-NE1-39) and 28 common genotypes (IT-CNNE1-28). Among these, 36 individuals were represented by a single genotype, while the remaining samples were grouped into genotypes shared by multiple individuals. Neutrality tests showed that Tajima’s D was negative but statistically insignificant, while Fu and Li’s D and F values were negative and statistically significant (P < 0.05). Additionally, Fu’s Fs was negative (Table S9). The mismatch distribution analysis revealed a unimodal pattern (Fig. S4).

All regions (ALL)

When all study areas were considered as a single population, 104 genotypes were identified. Among these, 37 genotypes were specific to CN, 39 were specific to NE, and 28 were shared between the two regions. Of the samples, 99 individuals were represented by a single genotype, while the remaining individuals were grouped into genotypes shared by multiple samples. Neutrality tests indicated that Tajima’s D was negative but statistically insignificant, while Fu and Li’s D and F values were negative and statistically significant (P < 0.02), and Fu’s Fs was also negative (Table S9). The mismatch distribution analysis produced a unimodal-like pattern (Fig. S4). Additionally, the analysis of polymorphic regions was conducted both for the combined population and for the two regional populations, with detailed results presented in Table S10.

Genetic differentiation and population structure

To assess genetic variance and pairwise genetic differentiation (F_ST) at the ITS2 gene level, populations were grouped into 31 sampling sites-based and two regional populations. AMOVA was performed locus by locus. In the dataset with 31 sampling sites, genetic variation among populations accounted for 7.31%, while within-population variation was 92.69%, resulting in an F_ST value of 0.07309 (P < 0.01). In the dataset divided into two regions, genetic variation among populations was calculated at 2.97%, while within-population variation was 97.03%, with an F_ST value of 0.02966 (P < 0.001) (Table S11). The most genetically differentiated populations, with the highest F_ST value (0.3188, P < 0.001), were observed between L10 in CN and L22 in NE. The distance matrix and color plot based on F_ST values for the 31 populations are provided in the supplementary materials (Fig. S5 and Table S12).

STRUCTURE analyses were performed with simulations ranging from K = 2 to 10 groups (27 simulations total, with three runs per group). The most appropriate group number was determined to be K = 2 using pophelper, with K = 3 (the second highest delta K) also considered for evaluation (Fig. S6). Results were assessed based on both two and three ancestral groups. When the simulation was based on two ancestral groups using the ITS2 gene, no distinct population structuring was observed at either the sampling site or regional level. However, the simulation with three ancestral groups revealed weak population structuring based on study regions, though it was less pronounced compared to the results from the cox1 gene (Fig. S7). The correlation between genetic variation and geographical distance was evaluated using the Mantel test for 361 D. marginatus individuals characterized by the ITS2 gene with individual geographical coordinates. The test yielded an r-value of 0.04177, which was statistically significant (P < 0.05). Additionally, SAMOVA analysis was conducted on the same dataset using geographical coordinates grouped into 31 sampling sites. Simulations from K = 2 to 10 groups were performed, with K = 4 determined as the most appropriate group simulation based on the F_CT value, although all groups were evaluated. At K = 4, sampling sites were not geographically segregated strictly by region. Instead, one group consisted of samples from nine sampling sites in CN (Group 1: L1, L3, L6, L4, L8, L10, L11, L12, and L15). Another group formed from a single location in NE (Group 2: L16), and a third group also in NE (Group 3: L18). A fourth, shared group (Group 4) included samples from 20 sampling sites across both CN and NE (L2, L5, L7, L9, L13, L14, L15B, LA, LGr, LM, L17, L19, L20, L21, L22, L23, L24, L25, L26, and L27) (Fig. S8).

Genotype distribution, Genotype network, and Pairwise-Distance

The 104 characterized D. marginatus genotypes based on the ITS2 gene were distributed as follows: 37 genotypes unique to CN, 39 unique to NE, and 28 shared between regions. Among the CN-specific genotypes, IT-CN7 was the most common, represented by three individuals, while IT-CN6, IT-CN9, and IT-CN16 were each represented by two individuals; the remaining genotypes were represented by a single individual. In NE, IT-NE20, IT-NE27, and IT-NE37 were each represented by two individuals, with all other genotypes represented by single individuals. For the 28 shared genotypes, IT-CNNE5 was the most abundant, represented by 56 individuals across most sampling sites, except L3, L13, L16, and L24. This was followed by IT-CNNE1 (30 individuals) and IT-CNNE2 (23 individuals) in various sampling sites. Other notable shared genotypes, such as IT-CNNE8 (20 individuals), IT-CNNE13 (19 individuals), and IT-CNNE9 (16 individuals), displayed varying degrees of distribution across CN and NE. The remaining shared genotypes were represented by fewer individuals, typically across limited sampling sites. Detailed information on the genotypes and their distributions across the different sampling sites is presented in Table S13 and Fig. S9.

A pairwise-distance analysis was performed using the 104 identified genotypes. The most genetically distant genotypes were IT-CN28 from CN and IT-NE29 from NE. IT-CN28 was represented by a single individual (Dm-13910) from location L13, while IT-NE29 was represented by a single individual (Dm-2421) from location L24. Intraspecific genetic variation ranged from 0.09% to 1.55%, with partial pairwise differentiation observed among region-specific genotypes. Notably, genetic differentiation within NE genotypes was slightly more pronounced compared to CN. However, this differentiation was less distinct than that observed in the cox1 gene (Table S14 and Fig. S9).

The 104 characterized genotypes were analyzed for genetic similarity and specificity using BLAST. The analysis confirmed that none of the genotypes had identical matches in existing records (based on full or near-full length comparisons), establishing that all 104 genotypes were unique. Among the most similar sequences, with homologies ranging from 98.36% to 99.91%, were D. marginatus records from Romania (FN296269, FN296278, FN296273, FN296275), Iran (GQ144707), and Germany (S83081), as well as sequences of closely related species, such as D. niveus from Iran (GQ144706) and D. silvarum from China (JQ737110) (Table S15).

Population structure and demographic history based on the cox1 + ITS2 concatenated dataset

The cox1 + ITS2 concatenated dataset (total length: 1923 bp) of 361 characterized D. marginatus individuals was aligned using the MAFFT algorithm. Comparisons within the dataset revealed no insertions or deletions. Analyses were performed separately for each sampling site, as well as at the regional level (CN and NE) and for the entire population dataset.

Central Anatolia (CN)

In the CN, 153 genotypes were identified, consisting of 149 region-specific genotypes (CNC-CN1-149) and four common genotypes (CNC-CNNE1-4) shared between regions. Among the samples, 129 individuals were represented by a single genotype, while the remaining individuals were grouped into genotypes shared by multiple samples. Neutrality tests revealed that Tajima’s D and Fu and Li’s D and F values were negative and statistically significant (P < 0.05), while Fu’s Fs was also negative (Table S16). The mismatch distribution analysis showed a unimodal pattern (Fig. S10).

Northeast Anatolia (NE)

In the NE, 129 genotypes were identified, including 125 region-specific genotypes (CNC-NE1-125) and four common genotypes (CNC-CNNE1-4) shared between regions. Among the samples, 121 individuals were represented by a single genotype, while the remaining individuals were grouped into genotypes shared by multiple samples. Neutrality tests showed that Tajima’s D was negative but statistically insignificant, while Fu and Li’s D and F values were negative and statistically significant (P < 0.05), and Fu’s Fs was also negative (Table S16). The mismatch distribution analysis resulted in a unimodal graph (Fig. S10).

All regions (ALL)

When all study areas were considered as a single population, 278 genotypes were identified. Of these, 149 genotypes were specific to CN, 125 were specific to NE, and four genotypes were shared between the two regions. Among the samples, 250 individuals were represented by a single genotype, while the remaining individuals were grouped into genotypes that were shared by multiple samples. Neutrality tests indicated that Tajima’s D and Fu and Li’s D and F values were negative and statistically significant (P < 0.02), while Fu’s Fs was negative (Table S16). The mismatch distribution analysis showed a unimodal pattern (Fig. S10). Additionally, the results from the analysis of polymorphic regions, considering both the combined population and the two regional populations, are presented in Table S17.

Population structure

To assess genetic variance and pairwise genetic differentiation (F_ST) within and between populations using the concatenated dataset, populations were categorized into 31 sampling sites-based and two regional populations. AMOVA was performed locus by locus. In the dataset with 31 populations, the genetic variation among populations was 10.51%, while within-population variation was 89.49%, resulting in an F_ST value of 0.10511 (P < 0.001). In the dataset based on two regions, genetic variation among populations was 9.4%, and within-population variation was 90.6%, with an F_ST value of 0.09396 (P < 0.001) (Table S18). The most differentiated population, with the highest F_ST value (0.28249, P < 0.001), was found between L14 in CN and L16 in NE. The F_ST distance matrix and color plot based on the 31 populations are provided in the supplementary files (Fig. S11 and Table S19).

STRUCTURE analyses were performed with simulations ranging from K = 2 to 10 groups (27 total simulations, three for each group). The optimal number of groups was determined to be K = 2 using pophelper, with K = 3 also considered due to the second highest delta K value (Fig. S12). The results were evaluated for both two and three ancestral groups. The simulation based on K = 2 showed less pronounced population structuring at both the sampling site and regional levels, based on the concatenated dataset. However, when K = 3 was used, regional population structuring became evident, though weaker than that observed for the cox1 gene (Fig. S13). To assess the correlation between genetic variation and geographical distance, a Mantel test was conducted on the dataset of 361 D. marginatus individuals, using individual geographical coordinates. The test yielded a statistically significant r-value of 0.1863 (P < 0.05). Additionally, SAMOVA analysis was performed with geographical coordinates based on 31 sampling sites. Simulations from K = 2 to 10 were carried out, with K = 5 determined as the most appropriate based on the F_CT value. This simulation revealed that, at K = 2, CN and NE were separated into two distinct groups. As the number of groups increased, separations were observed within the NE populations, while CN populations remained as a single group. From K = 6 onward, locations L9 and L14 were the first to separate from CN. At K = 5, the CN populations formed a single group, while the NE populations were divided into four groups: Group 1 (L17), Group 2 (L18), Group 3 (L19, L20, L21, L22, L23, L24, L25, L26, L27), and Group 4 (L16) (Fig. S14).

Genotype distribution, genotype network, and pairwise-distance

A total of 278 genotypes were identified in the characterized D. marginatus individuals based on the concatenated dataset. Among these, 149 genotypes were specific to CN, 125 genotypes were specific to NE, and four genotypes were shared between both regions. The genotypes from CN included the most common CNC-CN10, which was represented by nine individuals across multiple sampling sites (L1, L5, L8, L9, L10, L12, L15B, LGr), and CNC-CN14, which was also represented by nine individuals but across different locations (L2, L3, L4, L5, L6, L11, L12, L15). Other notable genotypes in this region were CNC-CN13, represented by eight individuals, and CNC-CN22, represented by seven individuals, among others. The remaining genotypes from CN were represented by fewer individuals, with some found in specific locations and others by just a single individual. In NE, the genotypes were similar, with the most common being CNC-NE9, CNC-NE38, and CNC-NE53, each represented by two individuals in specific sampling sites. The majority of the genotypes in NE were represented by one individual. The shared genotypes between CN and NE included CNC-CNNE2, represented by six individuals in multiple sampling sites, and CNC-CNNE1, which was found in three individuals across L3, L10, and L16, while other common genotypes showed variations across both regions. Detailed information on these genotypes and their distribution across the different sampling sites can be found in Table S20.

Pairwise-distance analysis was performed using the 278 genotypes identified in this study. Among these genotypes, the most genetically distant were CNC-NE93, CNC-NE94, CNC-NE15, and CNC-NE16, all of which were found in NE. Specifically, the CNC-NE93 and CNC-NE94 genotypes were represented by a single individual each in sampling site L24 (Dm-2411 and Dm-2421), while CNC-NE15 and CNC-NE16 genotypes were represented by single individuals in sampling site L17 (Dm-1714 and Dm-1721). The intraspecific genetic variation ranged from 0.05% to 1.20%, and pairwise differentiation among region-specific genotypes was observed, particularly in NE, where the differentiation within genotypes was more pronounced. These findings highlight regional genetic variation, and further details on the pairwise distance analysis can be found in Table S21 and Fig. S15.

Phylogenetic relationship

Cox1 phylogeny

A dataset including 131 D. marginatus haplotypes based on the cox1 gene was created, along with specimens recorded as D. marginatus, D. niveus, D. raskemensis, D. nuttalli, and D. silvarum in the subgenus Serdjukovia from the GenBank database. Additionally, three D. reticulatus specimens (MT478096, OM142141, and OQ947121), classified under a different subgenus of Dermacentor, were used as outgroups. After performing reliability analysis to remove duplicate sequences and ensure minimum length, the dataset was edited and aligned. The final dataset consisted of sequence data for a total of 281 samples (comprising 131 haplotypes, 147 GenBank records, and 3 outgroups) (Table S22).

The constructed phylogenetic tree revealed three main clades: D. marginatus, D. raskemensis, and D. silvarum/nuttalli, all supported by high posterior probabilities (Fig. 4). Within the D. marginatus clade, at least five haplogroups and three well-supported branches (posterior = 0.82–1) were identified, including 201 D. marginatus and seven D. niveus haplotypes. Haplogroup Dm1 was split into two subclades, with subclade Dm1a consisting of 44 haplotypes from this study (40 from NE, two from CN, and two common haplotypes) along with 23 D. marginatus haplotypes from Europe, Asia, and North Africa (Romania, Hungary, Croatia, Slovakia, France, Italy, Germany, China, Kazakhstan, and Tunisia). Subclade Dm1b included 86 haplotypes from this study (65 from CN, 19 from NE, and two common haplotypes) and 12 D. marginatus haplotypes from Europe, Asia, and North Africa (Romania, Croatia, Portugal, Türkiye, Iran, and Tunisia). The phylogenetic tree further revealed that haplogroup Dm2 consisted of a single haplotype from this study (CX-NE53) and 10 D. marginatus haplotypes from China and Kazakhstan. Haplogroup Dm3 included 14 D. marginatus haplotypes from Kazakhstan, China, and Russia. Haplogroup Dm4 comprised seven D. marginatus and seven D. niveus haplotypes from China and Kazakhstan, while haplogroup Dm5 consisted of two D. marginatus haplotypes from China and Kazakhstan. In addition, the D. marginatus main clade contains three separate branches supported by high posterior probabilities (posterior = 0.82–1) and including one haplotype each from Kazakhstan, Pakistan, and Iran. The D. raskemensis clade, with high posterior probability (posterior = 0.98), was divided into two subclades: one derived from a reference sequence of D. raskemensis from Türkiye (MT308586) and the other from a D. niveus sequence from Iran (MK863423). Finally, the D. silvarum/nuttalli main clade, which was highly supported (posterior = 0.98–1) and comprised 68 specimens, was divided into at least five distinct haplogroups. Haplogroup Dns1 included 23 sequences, comprising 16 D. nuttalli and seven D. silvarum haplotypes from China. Haplogroup Dns2 consisted of 29 D. nuttalli and six D. silvarum haplotypes from China, Mongolia, and Russia, while haplogroup Dns3 contained four D. nuttalli and two D. silvarum haplotypes from China. Haplogroup Dns4 comprised two D. nuttalli haplotypes from China, and haplogroup Dns5 contained two D. nuttalli haplotypes from China (Fig. 4 and Table S22).

Upon examining the phylogenetic positions of the haplotypes characterized in this study, the phylogenetic tree revealed that all haplotypes, except for CX-NE53, cluster into two subclades of haplogroup Dm1 alongside D. marginatus haplotypes from Europe, Asia, and North Africa. The haplotype CX-NE53, obtained from sampling site L24 in NE, clustered into haplogroup Dm2 alongside D. marginatus haplotypes from China and Kazakhstan. This haplogroup is monophyletic with haplogroup Dm3, another Central Asia-related haplogroup, supported by a maximum posterior probability (Fig. 4). The apparent separation of CX-NE53 from the other major haplogroup was further confirmed in a phylogenetic tree that included only the haplotypes from this study (Fig. S16a). Moreover, within the D. marginatus main clade, D. niveus records were found to cluster into a single haplogroup (Dm4), though these did not form a distinct group with strong posterior probability support. Notably, a single additional D. niveus record clustered with the D. raskemensis main clade, being monophyletic with the reference D. raskemensis specimen at maximum posterior probability. The phylogenetic relationship of two other closely related species in the subgenus Serdjukovia, D. nuttalli and D. silvarum, showed that the D. silvarum/nuttalli main clade divides into two clades: D. nuttalli only and D. nuttalli + D. silvarum, with strong posterior probability support. Further analysis revealed that the mixed group further splits into four clades with very high posterior probabilities (posterior = 0.98–1), where one clade consists solely of D. nuttalli individuals, while the remaining three contain records for both species (Fig. 4).

ITS2 phylogeny

The ITS2 dataset was constructed to include 104 D. marginatus genotypes characterized in this study, along with ITS2 sequence data from closely related taxa within the subgenus Serdjukovia (D. marginatus, D. niveus, D. raskemensis, D. nuttalli, and D. silvarum) retrieved from GenBank. To root the phylogenetic tree, ITS2 sequences from D. reticulatus (OM142152, OR428530, and S83080) were incorporated as outgroup sequences. The dataset was curated to remove duplicate sequences and align minimum lengths to match the study’s genotypes. After reliability checks, the final dataset consisted of 268 sequences (104 genotypes from this study, 161 sequences from GenBank, and 3 outgroup sequences). Bayesian inference phylogeny was conducted using the GTR + Γ model, identified as the optimal nucleotide substitution model for this dataset. The resulting phylogenetic tree was rooted using the designated outgroup sequences for proper evolutionary context. Details of the sequences used here are available in Table S23.

The phylogenetic analysis, after collapsing branches with low posterior probabilities, identified three main clades: D. marginatus, D. raskemensis, and D. silvarum/nuttalli. These were supported by high posterior probabilities, alongside one additional D. silvarum/niveus clade with low posterior support (Fig. S17). Within the D. marginatus main clade, two genogroups (DmIT1 and DmIT2) and one branch were identified, all supported by maximum posterior probabilities. This clade includes a total of 132 sequences: 128 D. marginatus, three D. niveus, and one D. silvarum. The largest genogroup, DmIT1, comprises 104 genotypes from this study, 20 D. marginatus genotypes from Iran, China, Romania, and Germany, two D. niveus genotypes from Iran and China, and one D. silvarum genotype from China. Additionally, a single D. marginatus genotype from China (D. marginatus isolate XJ058, KC203417) formed a distinct branch, monophyletic with DmIT1. Genogroup DmIT2, including three D. marginatus genotypes from China, was connected to these branches, forming part of the broader D. marginatus clade (Fig. S17). The D. raskemensis main clade, monophyletic to D. marginatus, is divided into two subclades with a moderate posterior probability (Posterior = 0.74). One subclade includes a reference sequence from Türkiye (D. raskemensis isolate IT-D1109, PP618825), while the other (DnmIT) comprises records from Iran, including three D. niveus and one D. marginatus (Fig. S17). A separate genogroup (DsnIT), containing four D. silvarum records from China and one D. niveus record from Iran, was externally linked to the D. marginatus and D. raskemensis clades with low posterior probability (posterior = 0.2). The phylogenetic analysis of the D. silvarum/nuttalli clade revealed a complex structure with 122 sequences organized into four genogroups and one distinct branch. The first genogroup, DnsIT1, is the largest, comprising 79 genotypes primarily identified as D. nuttalli, alongside a single D. silvarum sequence. These records were distributed across China, Mongolia, and Russia. The second genogroup, DnsIT2, included 24 D. nuttalli and four D. silvarum genotypes, also originating from China and Mongolia. The third genogroup, DnsIT3, was exclusively composed of nine D. nuttalli genotypes from Mongolia. The fourth genogroup, DnsIT4, contained six D. marginatus genotypes from China, while a separate clade consisted of a single D. nuttalli genotype (D. nuttalli isolate NM68, MW477873) from Mongolia, forming a distinct lineage (Fig. S17). Detailed information on sequences and genogroups can be found in Table S23.

The phylogenetic analysis of genotypes characterized in this study revealed that all genotypes clustered within the D. marginatus main clade, specifically into the genogroup DmIT1, which is the largest group in this clade and includes D. marginatus, D. niveus, and D. silvarum records from European and Asian specimens (Fig. S17). Within this genogroup, the genotypes were further divided into two distinct subgroups supported by maximum posterior probability. One subgroup consisted of 11 genotypes (four from NE, three from CN, and four shared genotypes), while the second subgroup contained 93 genotypes (34 from CN, 35 from NE, and 24 shared genotypes) (Fig. S16.b). Interestingly, D. marginatus records in the phylogenetic tree displayed notable paraphyly. The genogroup DnsIT4, comprising six D. marginatus records from China, was found to be paraphyletic with the D. marginatus main clade. Instead, this group was monophyletic with a D. nuttalli genogroup (DnsIT3) supported by a high posterior probability (posterior = 0.85), clustering into the D. silvarum/nuttalli main clade. Furthermore, D. niveus records did not form a distinct and separate clade with sufficient posterior probability support. Instead, D. niveus sequences were found in mixed clades alongside D. marginatus and D. silvarum records (Fig. S17).