Population Genetics of the Medicinal Leech Whitmania pigra Based on Mitochondrial Genes

1. Introduction

Leeches are annelids belonging to the class Hirudinea [1], a group of invertebrates widely distributed in diverse freshwater environments worldwide (e.g., ponds, rivers, and marshes), with a few species also inhabiting marine or moist terrestrial ecosystems. Globally, about 680 leech species have been recorded, and they play important ecological roles [2], acting as intermediate predators or parasites in material cycling and energy flow, while also possessing remarkable medical value. Leech saliva contains multiple potent bioactive components, among which hirudin is the most prominent [3,4]. This substance was first extracted by Haycraft in 1884 [5] and later named by Jacoby in 1904 [5]. As the most potent natural specific thrombin inhibitor known to date [6], hirudin can efficiently and precisely inhibit thrombin activity [7], preventing the conversion of fibrinogen to fibrin [8], thus exhibiting strong anticoagulant effects and being regarded as a representative natural anticoagulant [9]. Beyond its anticoagulant function, hirudin has multiple pharmacological potentials. Liu et al. [10] found that it exerts protective effects against diabetic nephropathy, nephrotic syndrome, and renal interstitial fibrosis in chronic kidney disease. Moreover, other leech extracts have been used for topical treatment of knee osteoarthritis, showing anti-inflammatory regulatory functions [11,12], and have demonstrated certain capacities to inhibit cancer cell proliferation and metastasis [13]. They also show application prospects in anti-fibrosis and regulation of metabolic disorders [14], treatment of herpes zoster [15], and promotion of skin wound healing [16]. In China, more than 90 leech species have been recorded, exhibiting a distinct south-to-north distribution gradient. They are mainly concentrated in the middle and lower reaches of the Yangtze River, Southwest China, and South China, where water systems are well-developed and climates are humid. In contrast, Northern populations are primarily found in the wetlands of Northeast China and localized areas of the Yellow River basin [17]. Among them, medicinal leeches have been extensively studied, such as Whitmania pigra and Hirudo nipponia, which are endemic or common species with significant economic and medicinal value [18,19,20,21].

W. pigra is one of the most widely distributed and medicinally valuable leech species in China, commonly found in freshwater bodies across various provinces, as well as in Russia, Japan, and other countries [22,23]. It has a moderate body size, typical suckers, well-developed jaws and teeth, and a complete digestive system [24,25,26]. It feeds on the body fluids of mollusks and prefers warm, shallow waters with stone and algae attachment [27]. More importantly, this species has been included in the Chinese Pharmacopoeia [28] and is currently the most widely used medicinal leech variety in clinical practice in China, with advantages such as accessible resources, strong ecological adaptability, and clear efficacy. It has thus become the core subject of research and utilization of leech medicinal materials, holding irreplaceable scientific and industrial importance [29,30,31,32,33].

Existing basic biological studies on this species are mostly focused on macroscopic aspects such as taxonomy, feeding habits, distribution, and behavioral traits, and remain relatively fragmented [34]. Taxonomically, the family placement of the genus Whitmania and the delimitation of closely related species are still controversial; traditional morphological classification cannot accurately distinguish sibling species, and the taxonomic system is not yet unified [35,36]. Feeding studies are mostly based on artificial breeding conditions, clarifying the effects of environmental factors on feeding and growth [37,38,39], but systematic investigations of long-term feeding adaptation mechanisms in natural populations are lacking. Distribution and behavioral studies have only provided a preliminary outline of macro-distribution patterns and basic activity rhythms; limited molecular marker work has been restricted to local populations, failing to elucidate population differentiation and behavioral adaptive evolution on a broader scale [40,41,42], and mechanistic studies at the population level are clearly insufficient. Conducting population genetics research holds significant theoretical and applied value. In application, it can clarify taxonomic and germplasm differences, providing a basis for medicinal material identification, quality standards, and medication safety; assess wild diversity to support conservation, guide selective breeding, and alleviate germplasm bottlenecks. In theory, it can explain population differentiation and adaptive mechanisms, analyze the genetic regulation of medicinal traits, and fill research gaps. In summary, this study is key to supplementing fundamental knowledge and promoting resource conservation and industrial development.

Mitochondrial DNA (mtDNA), as a molecular marker, is widely used in population genetics and phylogeography due to its unique genetic characteristics: strict maternal inheritance with no recombination during transmission, providing a clear and traceable genetic lineage for tracking maternal ancestry and population history [43]; a relatively high nucleotide substitution rate, accumulating abundant variable sites within species and among closely related species [44], sensitively reflecting population genetic diversity and differentiation patterns; high copy number and relatively conserved structure, allowing easy design of universal primers and high PCR amplification success rates, enabling stable information retrieval even from degraded samples [45,46]. The COI gene has been recognized as the standard DNA barcode for animals, and the Cytb gene is also a classic marker in population genetics studies of vertebrates and invertebrates [47,48,49], both serving dual functions of species identification and genetic diversity analysis, achieving “one marker for multiple uses”. Müller et al. [50] pointed out that mitochondrial protein-coding genes (including COI, Cytb, etc.) perform well in phylogenetic reconstruction of non-primate invertebrates; Tong et al. [51] also emphasized the potential application of such genes in exploring bioactive components in the W. pigra genome. However, mtDNA markers also have limitations, such as an effective population size only one-quarter that of nuclear genes, maternal inheritance failing to reflect male-mediated gene flow, and potential interference from nuclear mitochondrial pseudogenes (NUMTs) [52]. In practical applications, marker strategies should be carefully chosen and, when necessary, combined with nuclear markers for cross-validation [53].

Previous molecular population genetics studies on W. pigra, although using various markers, still have notable limitations. First, early studies relied on dominant markers with poor reproducibility and limited information content. Liu et al. [54] used ISSR and SRAP markers to analyze genetic diversity and detected a high proportion of polymorphic loci, but these dominant markers cannot distinguish homozygous from heterozygous states, making accurate estimation of genetic parameters difficult; amplification results are easily affected by PCR conditions, resulting in poor experimental reproducibility and difficult integration of data across laboratories. Second, geographic coverage was narrow, failing to reveal a national-scale differentiation pattern. Jiang et al. [55] used RAPD markers to analyze five populations from Xinghua (Jiangsu), Binxian (Heilongjiang), Dawu (Hubei), etc., and found significant regional differences, but sampling was still limited to parts of East and North China, without covering major southern distribution areas such as Sichuan, Hunan, Jiangxi, and Anhui. The generality of conclusions was constrained, insufficient for supporting national resource conservation planning. To date, no study has systematically depicted the genome-wide genetic structure of W. pigra across China. Third, single-gene marker strategies provide limited phylogenetic information. Yue et al. [56] used only the Cytb single gene to analyze three wild populations; although they detected 47 haplotypes with high haplotype diversity (Hd = 0.9781), the limited variable sites of a single gene cannot deeply reveal complex evolutionary relationships and differentiation histories. Hao et al. [57] expanded the sample to five populations and 146 individuals using COI and Cytb dual-gene markers, but were still restricted to two mitochondrial genes, lacking a more comprehensive genomic perspective.

This study adopted a combined three-gene (COI, Cytb, ND1) marker strategy, providing richer variable sites and phylogenetic information compared to single- or double-gene studies. Sampling covered 13 geographic populations across China with 137 individuals, involving Heilongjiang, Shandong, Anhui, Jiangsu, Hubei, Hunan, Jiangxi, Henan, Sichuan, Tianjin, and other major distribution areas, with substantially expanded geographic coverage and sample size relative to previous work. Multiple analytical methods were comprehensively applied to systematically dissect the genetic structure, differentiation patterns, and demographic history of the national populations. Collectively, this study achieves a systematic advancement over existing studies in terms of marker number, geographic coverage, sample size, and analytical approaches. The results will provide a solid scientific basis for formulating germplasm resource conservation strategies, defining Evolutionarily Significant Units (ESUs), and guiding selective breeding of W. pigra.

2. Materials and Methods

2.1. Sample Collection and Sequencing

A total of 137 wild W. pigra individuals were collected from 13 localities across China (Figure 1, Table 1). 10–12 individuals were collected from each population. Anterior tissue was excised from each individual, and total genomic DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen). Qualified DNA samples were used to construct libraries with ~350 bp insert fragments using Illumina-compatible library preparation kits, and 150 bp paired-end whole-genome resequencing was performed on the BGISEQ sequencing platform. Raw data were quality-controlled using FASTP v0.20.0 [58] to obtain clean reads for each sample.

2.2. Extraction of Mitochondrial Gene Sequences

The clean reads obtained from resequencing were assembled de novo using MEGAHIT v1.2.9 [59] to generate contig sequence files. This study focused on three mitochondrial protein-coding genes (COI, Cytb, and ND1) for population genetic analyses of W. pigra. Using the published mitochondrial genome of W. pigra as reference, the COI, Cytb, and ND1 sequences were used as bait sequences to search for homologous sequences in the assembled contigs via BLAST v2.13.0 [60]. The retrieved homologous sequences were aligned using MEGA v11.0.13 [61], and the homologous regions were accurately extracted.

2.3. Genetic Diversity Analysis

The coding region sequences of the three target mitochondrial protein-coding genes were integrated into separate FASTA files (strictly restricted to CDS regions). Subsequently, multiple sequence alignment was performed using the “Align by MUSCLE” function in MEGA v11.0.13. Based on the alignment results, DnaSP v6 [62] was used to calculate genetic diversity parameters at the population level for each gene: number of variable sites (S), number of haplotypes (Hn), haplotype diversity (Hd), and nucleotide diversity (Pi). The three genes were concatenated, and the same parameters were calculated for the combined dataset across all populations. To evaluate genetic variation within each geographic population, DnaSP v6 was also used to calculate the above parameters for each individual gene and the concatenated genes within each population.

2.4. Haplotype Networks and Phylogenetic Analysis

Based on the sequence datasets of the three genes (COI, Cytb, and ND1), haplotype network construction and phylogenetic analyses were performed separately for each gene. First, haplotypes were extracted from the alignment results of each gene using DnaSP v6. Then, HapSolutely v0.2.2 [63] was used to construct gene-specific haplotype networks under the Fitch model. Subsequently, homologous sequences of the closely related species Whitmania laevis (GenBank accession: NC_023926.1) were downloaded as Outgroup. MEGA11 was used to construct NJ phylogenetic trees for haplotypes of each gene and for the combined haplotypes. All analyses were performed with 1,000 bootstrap replicates to assess node support. The resulting tree files were visualized and optimized using FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/Figtree/). The same methods were applied to the concatenated gene dataset for haplotype network construction and phylogenetic analysis.

2.5. Analysis of Population Genetic Structure

Bayesian clustering analysis was performed using STRUCTURE v2.3.4 [64] for the three individual mitochondrial genes and the concatenated genes. The number of clusters (K) was explored from K = 2 to K = 7. Key parameters were set as follows: burn-in length of 50,000 steps, followed by 50,000 MCMC replicates after burn-in, with 20 independent iterations for each K value. CLUMPAK (https://clumpak.tau.ac.il/results.html) was used to analyze the results online and calculate ΔK to determine the optimal population structure (K) and generate genetic structure plots. Pairwise F_ST values between populations were calculated using Arlequin v3.5 [65]. Permutation tests were performed with 1,000 replicates, and the significance level (α) was set at 0.05.

2.6. Historical Simulation of Population Dynamics

The Bayesian Skyline Plot (BSP) model implemented in BEAST v1.10.4 [66] was used to reconstruct the historical effective population size (Ne) changes of W. pigra. Because different genes have different mutation rates, making time unification difficult, the three genes were concatenated for analysis to better explore the population demographic history. The optimal nucleotide substitution model was selected using jModelTest v2.1.10 [67]. Analysis parameters were configured in BEAUti v1.10.4 [66]: substitution model: partitioned model; base frequencies: estimated; molecular clock model: uncorrelated lognormal relaxed clock; its prior was set as fixed. Given the lack of direct data on mitochondrial gene evolutionary rates in annelids, we referred to the conservative substitution rate of insect mitochondrial protein-coding genes, set at 1.35% per site per million years.

3. Results

3.1. Variations in Gene Sequence

We obtained sequences of the target genes COI, Cytb, and ND1 from 137 individuals across 13 populations. Among them, COI had the longest length (1,534 bp), and ND1 the shortest (874 bp). Alignment revealed extensive genetic variation, with a total of 441 variable sites identified, ranging from 104 to 182 variable sites and 53 to 78 haplotypes per gene. COI had the highest number of variable sites and haplotypes, while ND1 had the lowest, indicating that sequence length substantially affects variation. For nucleotide diversity, Cytb ranked highest, followed by COI, while the shortest gene ND1 remained the lowest (Table 2).

3.2. Genetic Diversity Within a Gene Population

Using DnaSP v6, genetic diversity parameters were calculated for the 13 geographic populations of W. pigra based on COI, Cytb, and ND1 sequences (Table 3), including the number of variable sites (S), haplotype diversity (Hd), and nucleotide diversity (Pi), with combined information after concatenation.

Among the individual genes, COI exhibited the highest level of genetic variation, with 182 variable sites in total. JSTZ and HBWH had the highest number of variable sites (41 each), while JXJA had none (S=0). Haplotype diversity (Hd) was highest in JSTZ (1.000) and zero in JXJA, indicating no variation detected in the COI gene for that population. Nucleotide diversity (Pi) was generally moderate to high, consistent with the distribution of variable sites. For Cytb, total variable sites were 155, with TJ having the most (38). The highest Hd was in JSTZ (0.978) and the lowest in JXJA (0.200). Nucleotide diversity for this gene remained moderately high in most populations, suggesting a relatively high evolutionary rate. For ND1, total variable sites were 104, the lowest among the three, indicating stronger conservation. HBWH had the most polymorphic sites (18). Hd was highest in TJ (0.939) and still lowest in JXJA (0). Nucleotide diversity followed a similar pattern, further supporting the lower variability of ND1.

After concatenating the three genes, the total number of variable sites was 86, lower than the sum of individual gene counts, possibly reflecting non-additive variation across genes. JSTZ still exhibited extremely high genetic diversity in the combined dataset (Hd=1.000), highlighting its role as a reservoir of comprehensive genetic variation. highlighting its comprehensive accumulation of genetic variation. JXJA had a combined Hd of only 0.200, significantly lower than other populations, indicating a relatively homogeneous genetic background. Overall, different mitochondrial genes showed marked differences in genetic variation, with COI being the most variable and ND1 the most conserved. JSTZ consistently showed high or the highest genetic diversity across all three genes and the combined dataset, suggesting it may be a genetic diversity center for W. pigra with important conservation significance. Conversely, JXJA consistently showed low diversity, possibly having experienced a genetic bottleneck or founder effect, and thus requires priority genetic conservation measures.

3.3. Haplotype Network Analysis

The haplotype networks constructed based on the sequences of the three mitochondrial genes COI, Cytb, and ND1 all exhibit similar topological structures, clearly divided into three major branches (A, B, and C), with a characteristic radial distribution from the central haplotype outward (Supplementary Figures S1-S3). Specifically, Branch A consists primarily of haplotypes derived from AHBB, HBZX, SDLY, SDZZ, and TJ, along with some haplotypes from AHCZ, HNXY, and SCYB; Branch B comprises mainly haplotypes from JXJA, HBWH, JSTZ, and HNYY, including some from AHCZ and SCYB; whereas Branch C demonstrates significant geographic specificity, with all its haplotypes originating exclusively from the HLJ population. The combined haplotype network for the three genes (Figure 2) yields identical results to those derived from individual genes.

3.4. Phylogenetic Analysis

Using the Polar Tree layout visualization feature in FigTree, the phylogenetic trees of the three genes were transformed from a simple branch diagram into a circular “radar plot” (Supplementary Figures S4-S6). The haplotype-based phylogenetic trees for all three genes each exhibit three distinct branches (A, B, and C): Branch A primarily comprises the haplotypes AHBB, HBZX, SDLY, SDZZ, AHCZ, and TJ; Branch B consists mainly of JXJA, HBWH, HNYY, JSTZ, and SCYB; notably, Branch C contains only the haplotype HLJ. When the three genes were combined for phylogenetic analysis, the results aligned with those obtained from single-gene analyses (Figure 3), identifying three clonal lineages (A, B, and C). This phylogenetic structure corresponds precisely to the findings from haplotype network analysis.

3.5. Genetic Structure and Genetic Differentiation

The results of the three-gene STRUCTURE analysis, calculated online using CLUMPAK, yielded a consistent optimal ΔK value of 3 for all analyses, unequivocally supporting the division of the 13 geographic populations of China W. pigra into three genetically distinct core groups. The most significant finding from the genetic structure maps of COI, Cytb, and ND1 (Supplementary Figures S7-S9) is that the HLJ population consistently exhibited a highly homogeneous and independent genetic composition across all analyses, with individuals belonging almost entirely to a single ancestral component across all clusterings ranging from K=2 to K=7, demonstrating genetic isolation from other populations. Beyond this independent HLJ unit, the genetic structure map further revealed the geographical distribution of the remaining populations: one comprises southern populations centered on JXJA, JSTZ, HNYY, and HBWH, exhibiting close genetic correlations; the other includes northern-central populations such as TJ, SDLY, SDZZ, AHBB, and HBZX. Notably, certain central populations (e.g., AHCZ and SCYB) accounted for proportions across all three ancestral components, suggesting these regions may represent transitional zones or historical hybridization areas with some degree of gene flow or secondary contact. Combining the three genes for population structure analysis (Figure 4) yielded results consistent with those from single-gene analyses.

This genetic classification is strongly supported by quantitative results from pairwise genetic divergence (F_ST) analyses between populations. The pairwise genetic divergence data for the three genes individually (Supplementary Tables S7–S9) and when combined (Table 4) demonstrate that the F_ST values between the HLJ population and all other populations are exceptionally high (0.902–0.977, p < 0.001), indicating a level of differentiation approaching that of a subspecies and statistically confirming its extreme genetic uniqueness. Furthermore, the F_ST values between southern populations and northern-central populations are generally > 0.5, suggesting significant genetic differentiation among these groups as well.

These results are highly consistent with the previous haplotype network (all HLJ individuals forming an independent branch without shared mutational steps with other haplotypes) and phylogenetic analyses (HLJ forming a highly supported monophyletic group), constituting a comprehensive chain of mutually reinforcing evidence. In summary, the population genetic structure of Chinese W. pigra consists of three clear evolutionary units: an isolated and unique northeastern group (HLJ), a southern group, and a northern-central group. This pattern strongly suggests that historical geographic isolation events (e.g., glacial refugia isolation, major riverine barriers) have played key roles in the species’ evolution, and the HLJ group, as the most genetically distinct unit, should be considered an independent Evolutionarily Significant Unit (ESU) and prioritized for conservation.

3.6. Analysis of Population Historical Dynamics

The Bayesian Skyline Plot (BSP) results (Figure 5) clearly outline the population expansion–contraction history of W. pigra since the Late Pleistocene, tightly coupled with Quaternary climatic cycles. This history can be divided into four consecutive phases: P1 (0.2 Mya, late Middle Pleistocene), the population size remained at a low level with gentle fluctuations, reflecting a baseline stable state under multiple glacial–interglacial cycles; P2 (0.2–0.07 Mya), the population entered an unprecedented expansion phase, with exponential growth peaking during the Last Interglacial (LIG), when warm and humid climates greatly expanded its suitable habitats; P3 (0.07–0.012 Mya), the population size declined sharply due to harsh cold and dry conditions that reduced habitats; P4 (<0.012 Mya, since the Holocene), a rapid decrease associated with gradual global cooling and intensified anthropogenic impacts. Collectively, the population dynamics of W. pigra are the historical product of natural climatic oscillations and recent human disturbance.

4. Discussion

In this study, based on concatenated sequences of COI, Cytb, and ND1 (3,554 bp) from 137 W. pigra individuals across 13 geographic populations in China, a total of 441 variable sites were detected. All populations exhibited the combined pattern of high haplotype diversity (Hd > 0.5) and low nucleotide diversity (Pi < 0.02), This characteristic is widely interpreted as a typical signal of rapid expansion and accumulated mutations following a historical genetic bottleneck. Haplotype networks, phylogenetic trees, and population structure analyses consistently supported the division of the 13 populations into three significantly differentiated genetic groups: the northeastern group (HLJ), the southern group (JXJA, JSTZ, HNYY, HBWH), and the northern-central group (TJ, SDLY, SDZZ, AHBB, HBZX). Among them, the HLJ group exhibited extremely high genetic differentiation from all others (F_ST 0.902–0.977), suggesting it may represent an independent Evolutionarily Significant Unit (ESU). Bayesian Skyline Plot analysis reconstructed four phases of population demographic history, closely matching Late Pleistocene climatic fluctuations, indicating that glacial-interglacial cycles have profoundly affected population size changes in this species.

The pattern of high haplotype diversity (Hd) combined with low nucleotide diversity (Pi) revealed in this study is not unique among medicinal leeches. Lin et al. [68], based on 13 mitochondrial protein-coding genes from seven geographic populations of the Asian buffalo leech (Hirudinaria manillensis) in southern China, similarly found that all populations exhibited high haplotype diversity (> 0.5) and low nucleotide diversity (< 0.005), and accordingly inferred that the species had experienced historical bottlenecks followed by rapid expansion and mutation accumulation. Notably, that study further revealed that haplotype diversity increased logarithmically with gene length (R² = 0.858, p < 0.001), while nucleotide diversity showed a nearly perfect alternating low–high pattern (Z = 2.938, p = 0.003), suggesting that the interaction between gene length and evolutionary rate must be fully considered when interpreting mitochondrial genetic diversity parameters, and also providing quantitative support for the necessity of the three–gene combined strategy adopted in our study. Trontelj and Utevsky [69] showed, in their phylogeographic study of three medicinal leech species of the genus Hirudo, that despite wide distributions, mitochondrial diversity was generally low, implying Pleistocene bottlenecks and subsequent rapid recolonization from different glacial refugia following selective sweeps. Popa et al. [70], using mitochondrial COI and 12S markers on 133 samples of Hirudo verbana from Romania, further confirmed that this species is currently in a phase of population expansion, with wetland coverage and altitude being the primary ecological variables affecting its distribution. These comparative data indicate that medicinal leeches generally share similar bottleneck-expansion histories, but the phylogeographic fragmentation in W. pigra in this study is far more pronounced, with F_ST values between HLJ and other populations (0.902–0.977) far exceeding the inter–population differentiation levels observed in Hirudo species. This difference can be attributed to niche differentiation: as a non–blood–eeding leech, W. pigra lacks the passive dispersal pathways relying on host movement; its limited active dispersal capacity allows geographic isolation to accumulate deeper phylogeographic sorting.

The formation of the three–group genetic structure can be reasonably explained within the framework of geographic barriers, niche differentiation, and historical climatic events. The HLJ population consistently appeared as a highly homogeneous and independent unit in all STRUCTURE clustering schemes (K=2 to K=7), with its differentiation level reaching or exceeding that observed among many leech species, strongly suggesting that it is an independent ESU. The chromosome-level genome of W. pigra published by Liu et al. [71] identified 20 antithrombotic gene families; comparative genomic studies by Zhao et al. [72] on three non-blood-feeding Whitmania species showed that although all three have lost blood–feeding habits, they exhibit significant differences in the retention and evolution of antithrombotic–related genes, implying that different geographic populations may have already accumulated considerable genomic divergence, and the extreme differentiation of the HLJ group may be approaching the threshold of subspecific divergence. The differentiation between the southern group and the northern-central group highlights the barrier effect of the Qinling-Huaihe line as a key biogeographic divide in East Asia. Gao et al. [41], based on SNP markers, similarly found that five populations could be traced back to three population sources, qualitatively echoing the three major groups identified in our study. Hao et al. [57] (note: this should be ref. 57 in the original; the citation in the text appears to have a minor numbering error) based on COI and Cytb dual–gene analysis, showed that the Beijing population had the highest genetic diversity, while the Tianjin population had the lowest, with phylogenetic trees showing mixed distribution and no significant population expansion signal detected, largely consistent with our finding that the northern–central group (including TJ) had overall low diversity and lacked recent expansion signals, suggesting that this regional population may have maintained a small effective population size over a long period or experienced genetic homogenization due to aquaculture introductions. Oceguera-Figueroa [73] demonstrated that geographic barriers profoundly shaped speciation patterns in the blood-feeding leech genus Haementeria; the Qinling–Huaihe line, as a dual transition zone of climate and wetland ecosystem types, has a barrier effect highly parallel to those cases.

The four-phase population dynamics revealed by Bayesian Skyline Plot analysis precisely correspond to climatic cycles since the Late Pleistocene. The Last Glacial Maximum (LGM, approximately 26.5–19 ka) led to a substantial contraction of freshwater wetland habitats, causing a sharp reduction in the effective population size (Ne) of W. pigra and a genetic bottleneck; post-glacial warming and wetland expansion allowed remnant populations to radiate outward from multiple differentiated refugia. This process elucidates the mechanism underlying the high Hd-low Pi pattern: bottleneck events eliminate many ancient haplotypes, while rapid expansion allows a few surviving lineages to radiate numerous closely related haplotypes over a short period, but insufficient time to accumulate high levels of nucleotide divergence. Lin et al. [68] similarly revealed a five-phase demographic trajectory for Hirudinaria manillensis that closely matched climatic events, further confirming that glacial-interglacial cycles have profoundly affected the effective population sizes of leeches. The extreme genetic uniqueness of the HLJ population in our study suggests that it may have survived the LGM in a northern micro-refugium independent of southern refugia, and after deglaciation, it failed to undergo effective gene flow with southern lineages. The hermaphroditic, cross-fertilizing reproductive mode of W. pigra and its strict dependence on wetland habitats further intensified reproductive isolation and genetic drift among refugia, promoting pronounced phylogeographic divergence.

Taken together, our study clearly demonstrates that Chinese W. pigra populations exhibit a significant bottleneck-expansion genetic footprint, and the populations can be divided into three highly differentiated evolutionary clades—northeastern, southern, and northern-central—with the HLJ population as an independent ESU that should be prioritized as a conservation management unit. The Qinling–Huaihe geographic barrier and historical isolation in Northeast China are the main spatial drivers of lineage divergence, while Late Pleistocene glacial–interglacial cycles are the primary temporal drivers of population oscillations. This study not only fills the gap in population genetics research on this species based on combined mitochondrial multi-gene analysis, but also provides a new empirical case for the response paradigm of East Asian aquatic invertebrates to Late Pleistocene climate change, and has fundamental guiding value for genetic assessment, sustainable utilization, and selective breeding of medicinal leech resources. However, this study has certain limitations. First, analyses were based solely on mitochondrial protein-coding genes, which, as maternal markers, can only reflect maternal lineage history and cannot fully capture the overall population genetic characteristics of the species; future studies should incorporate nuclear markers such as microsatellites and SNPs to validate our conclusions from a biparental inheritance perspective. Second, although our sampling covers the main distribution areas of W. pigra in China, sampling density in key geographic transition zones (such as both sides of the Qinling–Huaihe line) remains insufficient, potentially affecting accurate resolution of gene flow patterns among populations; future studies should increase sampling density in transition zones and employ landscape genetic methods to quantify the barrier effects of geographic features on gene flow. Third, Bayesian Skyline Plot analysis based on mitochondrial gene mutation rates carries some uncertainty in time estimation; future work could combine paleoclimatic data and multiple independent molecular marker mutation rates for cross-validation. Fourth, this study did not involve morphological, ecological, or functional genomic investigations; future multidisciplinary approaches including morphometrics, ecological niche modeling, and transcriptomics could comprehensively elucidate the adaptive evolutionary mechanisms and conservation genetic needs of this species. Finally, given the special status of the HLJ population as an independent ESU, future research should focus on its current population status, habitat protection, and genetic resource management, and explore whether it meets the taxonomic criteria for species or subspecies revision, thereby providing a more solid scientific foundation for the conservation and sustainable utilization of W. pigra.

5. Conclusions

This study systematically elucidates the genetic structure of 13 geographic populations of W. pigra across China based on concatenated mitochondrial COI, Cytb, and ND1 gene sequences. All populations exhibit a typical pattern characterized by high haplotype diversity and low nucleotide diversity, which indicates rapid demographic expansion following historical population bottlenecks. Consensus results from haplotype network reconstruction, phylogenetic inference, and population structure analyses robustly support the delineation of three significantly differentiated clades, namely the Northeast, Southern, and Northern-Central groups. Notably, the Northeast (HLJ) population possesses extremely high genetic distinctiveness (F_ST>0.9), qualifying it as an independent Evolutionarily Significant Unit (ESU). Bayesian skyline plot (BSP) analyses reveal that demographic dynamics are tightly coupled with glacial-interglacial cycles during the Late Pleistocene, demonstrating that climatic fluctuations constitute the predominant historical driver of lineage differentiation. This research provides pivotal molecular evidence for the conservation of germplasm resources, formal definition of ESUs, and sustainable utilization of W. pigra. Future research shall integrate nuclear genomic markers to further dissect its adaptive evolutionary mechanisms and formulate targeted conservation strategies.