Field sampling was conducted as part of a survey of the malacofauna of the northern Tien Shan by the authors of this study (Ivan Nekhaev and Anel Ishayeva). The material examined was collected at two localities within the Lake Balkhash drainage basin in southeastern Kazakhstan, located approximately 200 km apart (Fig. 1).

At both localities, Radix plicatula was found to co-occur with the invasive species Physella acuta, and no other gastropod species were recorded.

The collected mollusks were examined using a Leica M125 stereomicroscope. Dissection of the specimens followed the methodology described by Kruglov (2005). The terminology of the reproductive system follows Jackiewicz (1993).

The collected specimens are deposited in the malacological collection of the Institute of Zoology (Almaty, Kazakhstan) under catalog numbers IZK-Mol-466, IZK-Mol-468, IZK-Mol-471, and IZK-Mol-478.

DNA was extracted from small pieces of foot tissue using the Qiagen DNeasy Blood and Tissue Kit (Qiagen, Netherlands), following the manufacturer’s protocol. A fragment of the mitochondrial gene cytochrome c oxidase subunit I (COI) was amplified using AllTaq Master Mix (Qiagen, Netherlands) with the primers LCO1490 and HCO2198 (BioBeagle, Russia), according to the protocol of Folmer et al. (1994).

The resulting PCR products were verified by agarose gel electrophoresis and sequenced in both directions using the Sanger method. Sequencing was performed by AllGene (Kazakhstan). All other laboratory procedures were carried out at the Institute of Zoology (Almaty, Kazakhstan).

The obtained sequences were deposited in GenBank; the accession numbers are available in Suppl. material 1: table S1.

For phylogenetic reconstruction, reference sequences from GenBank employed in the study by Saito et al. (2021), as well as nine newly generated sequences, were used. Sequences containing ambiguous nucleotide codes were excluded. The final dataset used for phylogenetic reconstruction comprised 65 sequences, including Galba truncatula (O. F. Müller, 1774) as the outgroup and other species of the genus Radix Montfort, 1810 (Suppl. material 1: table S2). For the construction of the Templeton, Crandall, and Sing (TCS) haplotype network, 253 sequences identified as Radix plicatula were used (Suppl. material 1: table S3).

Alignments and phylogenetic reconstructions were obtained using a previously described pipeline (Nekhaev 2023; Nekhaev et al. 2025b). Sequence alignment was performed using MAFFT. The optimal model of nucleotide substitution was selected in R and identified as TPM1u+G+I. The final alignment length was 469 nucleotides for phylogenetic tree reconstruction and 471 nucleotides for the TCS haplotype network.

The robustness of the resulting maximum likelihood (ML) tree topology was assessed using 100,000 bootstrap replicates. For phylogenetic reconstruction using Bayesian inference (BI), two independent runs were performed for 10 million generations, each consisting of one cold and three heated chains, with sampling every 1,000 generations. A burn-in of 2 million generations was applied, with parameter tuning performed every 100 generations.

Haplotype diversity (h), nucleotide diversity (π), and Tajima’s D were calculated in R using the same alignment as that used for the haplotype network analysis. Fu’s Fs neutrality test and mismatch distribution analysis, including the sum of squared deviations (SSD) and Harpending’s raggedness index, were performed using the Linux version of Arlequin v. 3.5 (Excoffier and Lischer 2010), executed in a Docker container.

The approximate mitochondrial lineage age was assessed using uncorrected pairwise COI distances among sequences belonging to the Tien Shan population. Pairwise p-distances were calculated, and only comparisons among Tien Shan sequences were retained for subsequent analyses. Mean and maximum within-lineage p-distances were then calculated from the upper triangle of the distance matrix.

To evaluate the sensitivity of lineage age estimates to uncertainty in mitochondrial substitution rates, divergence times were approximated under three alternative COI pairwise divergence scenarios: 1.0%, 1.5%, and 2.0% sequence divergence per million years. These rates were selected to encompass both slower and faster evolutionary scenarios around previously published estimates for invertebrates, including a reported mitochondrial divergence rate of approximately 1.57% per million years (Wilke et al. 2009).

The approximate lineage age was estimated by dividing the observed p-distances by the assumed divergence rate. Both mean and maximum within-lineage distances were converted into temporal estimates, with the maximum distance interpreted as an approximation of the mitochondrial coalescence depth of the local lineage.

This analysis was intended to provide an approximate temporal framework for evaluating whether the observed mitochondrial diversity is compatible with a recent introduction or longer-term persistence in the region, rather than to produce a formal molecular dating estimate. Such an approach provides only a rough approximation of lineage age because precise molecular dating is difficult when based on a small number of substitutions in a mitochondrial marker (Papadopoulou et al. 2010; Duchene et al. 2020).

In addition to the localities described in the present study, occurrence data for modeling were obtained from GBIF (available at https://doi.org/10.15468/dl.vvb67w), comprising 208 records of the species, and from the dataset published by Saito et al. (2021), which includes 235 occurrence records of Radix plicatula. These datasets were merged into a single occurrence database, after which duplicate records within a specified distance were removed. Several thinning radii were tested: 0.01 km (183 occurrences after thinning; Fig. 1), 0.5 km (174 occurrences), and 5 km (149 occurrences).

Two types of environmental predictors were used for modeling. The first consisted of bioclimatic variables from the WorldClim database (Fick and Hijmans 2017) at a spatial resolution of 30 arc-seconds (approximately 1 × 1 km). The second included land-cover variables derived from the Climate Data Record (CDR) and Interim CDR Sentinel-3 Land Cover dataset (v2.1.1) for 2022 (Bartholomé and Belward 2005). Prior to analysis, all environmental layers were cropped to the study region bounded by the coordinates 63°N, 25°E; 63°N, 167°E; 6°N, 25°E; and 6°N, 167°E.

A set of climatic variables corresponding to the physiological constraints of freshwater mollusks was used, including four variables describing the temperature regime—Annual Mean Temperature (Bio1), Mean Diurnal Range (Bio2), Temperature Annual Range (Bio7), and Mean Temperature of Wettest Quarter (Bio8)—and four describing the precipitation regime—Annual Precipitation (Bio12), Precipitation of Driest Month (Bio14), Precipitation Seasonality (Bio15), and Precipitation of Coldest Quarter (Bio19). This set of variables corresponds to that used by Saito et al. (2021), with the exception of Precipitation of Warmest Quarter (Bio18), which showed a high Spearman’s rank correlation coefficient (ρ = 0.82) within the study area.

Land-cover classes were grouped into seven categories: Wetland and Flooded Vegetation (codes 160, 170, and 180), Open Herbaceous Vegetation (codes 110 and 130), Deciduous Forests (codes 60, 61, 62, and 90), Coniferous Forests (codes 70, 71, and 72), Shrubs (codes 120, 121, and 122), Sparse Dry Cover (codes 150, 151, 152, 153, 200, 201, and 202), and Anthropogenically Modified Areas (codes 10, 11, 12, 20, 30, 40, and 190). Because these land-cover layers have a higher spatial resolution than the climatic layers, they were resampled to match the WorldClim grid. For each aggregated cell, the proportional cover of each land-cover group was calculated.

Collinearity among all selected layers was assessed using Spearman’s rank correlation coefficient and did not exceed 0.75 (Suppl. material 1: fig. S1). Species distribution models were built separately using either climatic variables alone or the climatic variable set combined with land-cover variables.

To reconstruct the potential climatic niche in previous geological epochs, the same set of bioclimatic variables from PaleoClim (Brown et al. 2018) was used, with the exception of Bio2 and Bio7, which are unavailable for some periods. Data were used for all 10 available intervals: late Holocene, Meghalayan (4.2–0.3 ka); mid-Holocene, Northgrippian (8.326–4.2 ka); early Holocene, Greenlandian (11.7–8.326 ka); Younger Dryas Stadial (12.9–11.7 ka); Bølling–Allerød (14.7–12.9 ka); Heinrich Stadial 1 (17.0–14.7 ka); Last Interglacial (ca. 130 ka); Marine Isotope Stage 19 (ca. 787 ka); mid-Pliocene warm period (3.205 Ma); and Pliocene M2 (ca. 3.3 Ma).

A separate model was trained at a resolution of 2.5 arc-minutes (approximately 4.5 × 4.5 km), the maximum resolution available in PaleoClim, without using land-cover predictors and excluding unavailable variables. Because the ranges of predictor values in the paleoclimatic data often did not match those used for model training, pixels for which the number of mismatched variables exceeded the median value across the entire study area were excluded during the projection of mollusk distributions onto paleoclimatic conditions. This approach reduces the estimated potential ecological niche of the species but avoids excessive extrapolation of the model to conditions outside the training range.

Species distribution modeling was performed following Nekhaev et al. (2026) using four algorithms: Generalized Linear Models (GLM), Generalized Boosted Models (GBM), Random Forests (RF), and Maximum Entropy (MaxEnt). The outputs of the individual algorithms were subsequently combined into an ensemble model using a weighted mean approach, with weights proportional to the predictive accuracy of each model (Araujo and New 2007).

Following the recommendations of Barbet-Massin et al. (2012), 10,000 pseudo-absence points were generated within a buffer of a specified radius around the occurrence records. Several buffer sizes were tested: 10, 50, and 100 km, as well as an option in which the entire study area was used as the buffer.

Model performance was evaluated using three complementary metrics: the area under the receiver operating characteristic curve (AUC), the true skill statistic (TSS), and Cohen’s κ (Kappa), extracted from BIOMOD2 model outputs.

All Supplementary materials, including code, tables, and figures, are available at https://doi.org/10.5281/zenodo.19398115.

Alignment construction and filtering for phylogenetic reconstruction and the TCS haplotype network were performed in R using the packages ape (https://github.com/emmanuelparadis/ape), ips (https://github.com/heibl/ips), and seqinr (https://github.com/lbbe-software/seqinr). Selection of the optimal model of nucleotide substitution and reconstruction of the ML tree were carried out using the packages ape and phangorn (https://github.com/KlausVigo/phangorn), whereas the TCS haplotype network was constructed using the packages pegas (https://github.com/emmanuelparadis/pegas) and igraph (https://github.com/igraph/rigraph). Bayesian phylogenetic reconstruction was performed in RevBayes v. 1.2.5 (Höhna et al. 2016). Scripts available at https://github.com/anisus3/Phylogeny were used with the necessary modifications. The final versions of the scripts are available in the Suppl. material 1 at https://doi.org/10.5281/zenodo.19398115.

Raster layers used for modeling were prepared in R using the packages terra (https://github.com/rspatial/terra), geodata (https://github.com/rspatial/geodata), and corrplot (https://github.com/taiyun/corrplot). Species distribution modeling was performed using the biomod2 package (https://biomodhub.github.io/biomod2/). The modeling workflow was based on previously developed scripts available on GitHub (https://github.com/anisus3/ecological_niche_modeling_v2), with minor modifications. The final versions of the scripts used in this study are provided in the Supplementary material.

Map preparation was carried out in QGIS v. 3.4 (https://qgis.org), using the Natural Earth Map Kit (https://www.naturalearthdata.com) as a basemap. Image editing (background removal), preparation of photographic plates, and final layout of computer-generated figures were performed in Pixelmator Pro v. 3.16 (https://www.pixelmator.com/pro/) and Sketch v. 77 (https://www.sketch.com).

All analyses were conducted on a computer equipped with an Apple Silicon M4 Max processor using the respective software versions listed above.

The collected specimens were relatively small, with the shell height of the largest individual reaching 8.8 mm (Fig. 2A). The shell height is approximately twice the shell width. The height of the spire varies slightly among the examined specimens (Fig. 2) and constitutes approximately 0.25–0.33 of the aperture height.

The length of the penis sheath (Fig. 2, ps) is approximately equal to the length of the praeputium (Fig. 2, p). The pyriform body (Fig. 2, pb) is elongated and trapezoidal in shape. The spermatheca (Fig. 2, s) is large and connected to the vagina (Fig. 2, v) by a duct (Fig. 2, ds) whose length is nearly twice the diameter of the spermatheca.

Overall, the examined specimens correspond well to the descriptions and illustrations of Radix plicatula presented in recent literature (Vinarski et al. 2020; Aksenova et al. 2024).

Overall, the phylogenetic reconstructions are broadly consistent with previous studies (Aksenova et al. 2018; Vinarski et al. 2020; Saito et al. 2021) but provide substantially lower node support. In particular, the phylogenetic relationships of Radix alticola (Izzatullaev, Kruglov & Starobogatov, 1983) and Radix euphratica (Mousson, 1874) remain unresolved. This is most likely attributable to the limited ability of the COI fragment to resolve deeper phylogenetic relationships within the genus.

Both ML and BI analyses revealed substantial genetic heterogeneity within Radix plicatula. Some sequences from Taiwan and eastern mainland China were recovered as sisters to the Pamir species Radix alticola, whereas several sequences from eastern China clustered with Radix euphratica, a species distributed in the Middle East and the Pontic region. Although some of these clades received low support, they did not fall within the main clade containing the majority of Radix plicatula sequences.

Most sequences formed a single clade that was not supported in either the BI or ML tree (Suppl. material 1: figs S2, S3), including sequences from eastern China near the type locality of the species. This clade also contained the sequences obtained from the northern Tien Shan.

Radix plicatula exhibits high haplotype diversity, with a total of 118 haplotypes identified among 235 analyzed sequences (Fig. 3). The nine sequences obtained from the northern Tien Shan belonged to six distinct haplotypes, resulting in relatively high haplotype diversity for a geographically localized population (h = 0.89). Although this value was somewhat lower than that observed across the previously known range of the species (h = 0.99; 112 haplotypes among 226 sequences), nucleotide diversity within the Tien Shan samples also remained comparatively high (π = 0.65%), whereas overall nucleotide diversity across the entire dataset reached π = 2.9%.

The Tien Shan haplotypes were separated from the most closely related haplotypes from China, Japan, and Korea by two or three substitutions (Fig. 3), whereas one of the newly identified haplotypes differed from its closest relative by five substitutions. The Tien Shan haplotypes form a star-like pattern typical of groups sharing a common ancestor. However, this structure also includes haplotypes from geographically distant regions. Tajima’s D calculated for the Tien Shan population was close to zero and was not statistically significant (D = 0.173, p = 0.863). However, Fu’s Fs was strongly negative and highly significant (Fs = –6.839, p = 0.00017), indicating an excess of rare haplotypes relative to neutral expectations. The mismatch distribution was also consistent with the sudden demographic expansion model (SSD = 0.0215, p = 0.592; Harpending’s raggedness index = 0.0725, p = 0.721).

Pairwise comparisons between the Tien Shan and non-Tien Shan samples yielded a mean p-distance of 2.8% (range: 0.4–6.2%), indicating that the Tien Shan population harbors multiple divergent haplotypes rather than a single narrowly derived lineage.

Approximate mitochondrial lineage age estimates for the Tien Shan population varied depending on the assumed COI divergence rate. Mean within-lineage divergence corresponded to approximately 0.32–0.64 Myr, whereas estimates based on the maximum observed within-lineage divergence ranged from approximately 0.65 to 1.27 Myr.

A total of 19 different model variants were constructed using different combinations of buffer size, occurrence thinning, and predictor sets—climatic variables only and climatic variables combined with land cover. According to the performance metrics of the individual algorithms, the model built with thinning at 0.01 km, a buffer equal to the entire study area, and climatic predictors only (Model 1) performed substantially better than models constructed with restricted buffers (Suppl. material 1: figs S4–S5, table S4). These differences remained consistent in the ensemble evaluations, in which the model with a buffer equal to the study area showed the best performance (AUC = 0.99, TSS = 0.93, κ = 0.79) (Suppl. material 1: fig. S4, table S4). The second-best model used both climatic and land-cover predictors with a 50 km buffer and thinning at 0.01 km (Model 2) (AUC = 0.98, TSS = 0.90, κ = 0.60).

In Model 2, the contribution of climatic variables was consistently higher than that of land-cover predictors across all modeling algorithms (GLM = 85%, GBM = 82%, RF = 95%, and MaxEnt = 58%). In both models, the different algorithms evaluated the contribution of predictors differently; therefore, it was not possible to identify a consistent set of key variables across methods (Tables S5, S6).

The projections of the two models differed substantially. The projection of Model 1 closely matched that published previously (Saito et al. 2021), in which the potential niche was largely confined to East and Southeast Asia—the region from which Radix plicatula was previously known—with only small areas of suitable habitat in the eastern Pontic region (Fig. 4, Model 1). In contrast, the projection of Model 2 indicated much higher suitability across the study region, including northern Central Asia and a large part of Europe encompassed within the modeling extent (Fig. 4, Model 2). Two additional projections, built for a model with a thinning radius of 5 km using both climatic and land-cover predictors and for a model using only climatic predictors, were generally consistent with the projection of Model 2 (Suppl. material 1: figs S6, S7).

To project the present-day distribution of Radix plicatula onto paleoclimate, model parameters that provided the highest performance among the climate-only models at a resolution of 30 arc-seconds were selected; this was the model with a 100 km buffer and thinning at 0.01 km (Suppl. material 1: table S4). This model was then recalculated at a resolution of 2.5 arc-minutes using six climatic variables. The performance of this model was AUC = 0.97, TSS = 0.83, and κ = 0.52.

A large proportion of pixels in all paleoclimatic projections contained at least one predictor variable with values outside the range used for model calibration. The median number of such variables per pixel was 1.8–1.9 across all time periods, indicating widespread extrapolation. To limit its influence, all pixels in which two or more variables exceeded the training range were excluded. This threshold was chosen as a balance between retaining spatial coverage and minimizing extrapolation.

All paleoclimatic projections differed from the present-day projection (Fig. 5, Suppl. material 1: fig S9–S18). The main difference in most cases was the presence of an almost continuous area of high habitat suitability for Radix plicatula north of the major Asian mountain systems, which is absent under present-day conditions.

The previously known range of Radix plicatula includes southern China, Vietnam, Laos, Thailand, and the southern part of the Russian Far East (Brandt 1974; Aksenova et al. 2018; Vinarski et al. 2020; Saito et al. 2021). The specimens from the Tien Shan examined in the present study were identified as Radix plicatula based on both morphological characters and genetic marker analysis. This represents the first reliable record of this species outside its previously known range in East and Southeast Asia, approximately 2,200 km from the nearest documented locality.

Previously, a record of conchologically similar specimens identified as Pseudosuccinea casta (I. Lea, 1841) was reported from Kazakhstan (Vinarski et al. 2015). That record was based on an empty shell found in northwestern Kazakhstan, approximately 1,500 km from the locality reported here. The shell of that specimen is similar to those of Radix plicatula described in the present study. However, additional investigations are required to confirm the identity of that specimen and to determine whether the species has established a stable population in the region.

Radix plicatula exhibits a highly diverse phylogeographic structure both within its main range and in the Tien Shan. Some strongly divergent populations from Taiwan and China likely represent separate species. Clarifying their taxonomic status requires a dedicated study.

A previous study on the distribution and dispersal of Radix plicatula suggested, based on modeling of its potential range, that it is constrained by the contemporary climate, as the model projection largely corresponded to the known distribution (Saito et al. 2021). The authors applied an approach without occurrence thinning, apart from duplicate removal, and without restricting the buffer used for background selection in the MaxEnt algorithm. This approach is broadly comparable to the model used in the present study with thinning at 0.01 km and pseudo-absence selection across the entire study area. Despite differences in the modeling framework—the use of an ensemble approach rather than a single algorithm and the incorporation of a larger number of occurrence records, including GBIF data and the records reported here—the resulting projections are nearly identical to those obtained in the previous study. Importantly, Central Asia, including the foothills of the northern Tien Shan where the records were obtained, does not fall within the predicted suitable area in either the previous study (Saito et al. 2021) or Model 1 of the present study (Fig. 4, Model 1).

Although this model showed the highest performance metrics, it may provide an inaccurate representation of the species’ potential ecological niche. Its main limitation is the use of a buffer that coincides with the full modeling extent. The buffer defines the area within which pseudo-absences, or background points in the case of MaxEnt, are generated for comparison with presence records. Buffer selection should account for the ecological and biogeographic characteristics of the species, for example, its dispersal ability, and is typically restricted to intermediate spatial scales rather than the entire study region (VanDerWal et al. 2009; Barbet-Massin et al. 2012). Otherwise, the model primarily learns to distinguish the region of current occurrence from all other areas based on predictor combinations, leading to artificially inflated performance metrics without improving biological realism (Lobo et al. 2008). Despite this limitation, such an approach is commonly used in modeling the distributions of freshwater mollusks (Stensgaard et al. 2013; Daniel et al. 2018; Yang et al. 2018; Vikhrev et al. 2024).

Another potential limitation affecting both models presented in Fig. 4 is the relatively small thinning radius. The use of closely spaced occurrence records may increase spatial autocorrelation in the data, leading to biased parameter estimates and overestimation of model accuracy (Dormann et al. 2007; Veloz 2009). However, in the case of Radix plicatula, changing the thinning radius only slightly affects the number of records used for modeling and therefore has a limited impact on the results; the projection for a model with a thinning radius of 5 km is provided in Suppl. material 1: fig. S5.

The model with a restricted buffer (Fig. 4, Model 2) is therefore considered not only to better reflect the observed distribution of the species but also to provide a more accurate representation of the present-day potential ecological niche of Radix plicatula.

The much smaller current range of Radix plicatula compared with that predicted by the model (Fig. 4, Model 2) may be explained by both limited dispersal ability and competition with ecologically similar and widespread species of the genera Radix and Ampullaceana. Notably, in the Caucasus–Pontic region and southern Eastern Europe, where the model predicts high suitability for Radix plicatula, a phylogenetically close species, Radix euphratica (Mousson, 1874), has been recorded (Vinarski et al. 2020). Similar cases, in which genetically divergent lineages occupy environmentally similar regions but do not overlap in their distributions, have also been documented in other members of the family (Von Oheimb et al. 2011; Nekhaev et al. 2025b).

The main question of the present study is whether the record of Radix plicatula represents the result of a recent introduction, for example, through human-mediated dispersal, or the discovery of a population that has persisted in the region for a long time.

The principal pathways for the introduction of invasive freshwater mollusks into Central Asia are accidental or deliberate introductions into large water bodies and reservoirs associated with fisheries, as well as unintentional releases from the aquarium trade (Babushkin et al. 2023; Nekhaev et al. 2025a). Through these pathways, several taxa of Southeast Asian origin have entered the region, including Corbicula fluminalis (O. F. Müller, 1774), Sinanodonta woodiana (Lea, 1834), and Sinanodonta lauta (von Martens, 1877), as well as the Ponto-Caspian Monodacna colorata (Eichwald, 1829) and the North American species Physella acuta (Draparnaud, 1805) (Izzatullaev 2019; Babushkin et al. 2023; Nekhaev et al. 2025a).

Although Physella acuta was the only other gastropod species found together with Radix plicatula at both newly discovered localities, this association probably does not reflect a specific ecological relationship. Physella acuta is currently one of the most widespread freshwater gastropods in southern Kazakhstan and commonly occurs in a broad range of natural and artificial water bodies together with numerous native mollusk species (Nekhaev et al. 2025a; Nekhaev and Ishayeva 2026). Instead, the low species richness observed at both sites is more likely related to the characteristics of the habitats themselves. The pasture pond is a small temporary water body sampled under late-autumn conditions, whereas the hydrothermal spring represents an artificial warm-water habitat. Both environments are expected to support relatively species-poor mollusk assemblages.

These pathways of introduction appear unlikely for Radix plicatula, as the species has not been recorded in large water bodies where deliberate introductions associated with fisheries have occurred, for example, in Lake Balkhash (Nekhaev et al. 2025a; Nekhaev and Ishayeva 2026). In addition, this species is not commonly kept in the aquarium trade.

Invasive or recently established populations of freshwater mollusks typically exhibit pronounced founder effects, characterized by the dominance of one or a few haplotypes and substantially reduced genetic diversity relative to native populations. Such patterns often allow identification of the source population and, in some cases, estimation of the timing of introduction (May et al. 2006; Pigneur et al. 2011; Fontanilla et al. 2014; Albrecht et al. 2025; Nekhaev et al. 2025b; Schächinger et al. 2025).

A clear example is provided by the introduced South American populations of Galba truncatula, in which only four COI haplotypes have been recorded despite evidence for at least two independent colonization events originating from genetically distant western European and eastern Eurasian source populations (Nekhaev et al. 2025b). The invasive populations retain only a small fraction of the diversity present across the native range. Haplotype and nucleotide diversity in the native populations reach h = 0.90 and π = 3%, respectively, whereas the introduced populations exhibit markedly reduced variation. Similarly, most introduced populations of Pseudosuccinea columella worldwide belong to a single mitochondrial haplotype and a single multilocus genotype (“worldwide invader”), indicating an extreme reduction in genetic diversity during invasion (Lounnas et al. 2017).

Another comparison may be provided by the Icelandic lineage of Ampullaceana balthica, whose age was estimated at approximately 5 kyr and is therefore inconsistent with a recent anthropogenic introduction (Bolotov et al. 2017). Despite its likely natural origin and greater age, this lineage still exhibits signatures of founder events, including relatively low genetic diversity (h = 0.85, π = 0.3%) and a pronounced star-like haplotype network. Although two evolutionary lineages were identified in Iceland, the diversity within each lineage was even lower. Thus, founder-effect signatures may persist for several thousand years following colonization.

The northern Tien Shan population of Radix plicatula shows some similarities to the Icelandic case. It also contains a partially star-like haplotype structure and exhibits moderate levels of genetic diversity (h = 0.89, π = 0.65%). This genetic pattern is further supported by a strongly negative Fu’s Fs and a mismatch distribution consistent with the sudden demographic expansion model. However, unlike most documented invasive populations, it contains six haplotypes among only nine sequenced individuals, including several relatively divergent haplotypes. Furthermore, the star-like component of the network includes haplotypes from geographically distant regions rather than forming an isolated regional cluster. This pattern may reflect either historical bidirectional exchange between Central and East Asian populations or incomplete sampling of genetic diversity across the main range of the species. Even if future studies demonstrate that some of the haplotypes detected in the Tien Shan are not unique to this region, the observed level of diversity appears more consistent with a long-established population than with a recent introduction followed by rapid local diversification.

Some invasive species, including mollusks, may exhibit the so-called “genetic invasion paradox,” a pattern in which genetic diversity within introduced populations is comparable to or even exceeds that observed in the native range (Roman and Darling 2007; Dlugosch and Parker 2008; Hufbauer 2008). In most documented cases, this pattern is associated with multiple introductions from different source populations, resulting in admixture and the formation of novel combinations of pre-existing genetic variants rather than entirely new evolutionary lineages. Such processes may increase adaptive potential and accelerate evolutionary responses in introduced populations. However, invasive populations demonstrating this pattern usually retain a clear genetic signal of contributions from geographically distant native lineages (Dlugosch and Parker 2008), which is not observed in the present dataset. The Tien Shan haplotypes of Radix plicatula instead form a compact group closely related to East Asian lineages and separated from them by only a few substitutions.

It is also important to note that the “genetic invasion paradox” has primarily been documented in widespread and highly invasive species, so-called “global invaders,” often associated with repeated large-scale anthropogenic introductions. Radix plicatula does not fit this pattern, as it is currently known as a comparatively localized East and Southeast Asian species lacking evidence of extensive invasive expansion across multiple continents or climatic regions. Therefore, although multiple introductions could theoretically explain the relatively high haplotype diversity observed in the northern Tien Shan populations, this scenario appears unlikely because all detected haplotypes remain closely related to one another and to East Asian lineages, whereas the overall genetic pattern is consistent with demographic expansion within a single evolutionary lineage.

Under all tested evolutionary scenarios, the observed mitochondrial diversity was therefore compatible with a Middle–Late Pleistocene time scale rather than with a very recent introduction. This interpretation is not contradicted by the demographic analyses, which indicate population expansion but do not provide an absolute time scale for this event. This estimate is consistent with divergence times reported for other freshwater mollusk taxa, in which East Asian or Tibetan and Central Asian populations diverged during the Early and Middle Pleistocene (Von Oheimb et al. 2011, 2013).

The principal Late Pleistocene dispersal route from Southeast and East Asia to Central Asia discussed in the recent literature involves a series of unconnected lakes and river basins across the Tibetan Plateau. Such a pathway has been proposed for several groups of freshwater organisms, including Radix auricularia, a species that is both phylogenetically and ecologically close to Radix plicatula (Von Oheimb et al. 2011; Clewing et al. 2016; Deng et al. 2021). At the same time, the fragmented hydrographic structure of the plateau has promoted isolation among freshwater populations, leading to independent evolutionary trajectories in different regions of Tibet (Von Oheimb et al. 2013). Nevertheless, this route is not considered the most likely explanation for the origin of the Central Asian populations of Radix plicatula. Despite intensive molecular studies of Lymnaeidae conducted across the Tibetan Plateau during the last two decades, the species has never been recorded from Tibet, except along its extreme southwestern margin (Von Oheimb et al. 2011; Clewing et al. 2016; Saito et al. 2021; Vinarski et al. 2022).

A more plausible scenario involves a Late Pleistocene, or possibly even Holocene, dispersal route north of the Tibetan Plateau and the associated mountain systems. This hypothesis is strongly supported by paleoclimatic projections of the potential distribution of Radix plicatula.

Interpretation of projections of the present-day species distribution onto paleoclimate is associated with several limitations. First, paleoclimatic models are based on a smaller number of climatic variables—six instead of the eight used for the present-day niche—and do not include predictors such as land cover, which improve model performance. Second, paleoclimatic reconstructions are themselves model-based and may contain uncertainty (Brown et al. 2018). Finally, the species’ ecological preferences may have changed over time, and the present-day niche does not necessarily fully reflect the conditions under which it existed in the past.

Taking these factors into account, the reconstructions obtained should be regarded as a hypothesis of potentially suitable areas rather than a direct reconstruction of the actual range. At the same time, the agreement between the modeling results and molecular phylogenetic data, which indicate likely dispersal of the species from East Asian populations during the Pleistocene, suggests that the existence of extensive areas north of the Tibetan Plateau and the Tien Shan suitable for Radix plicatula during the periods considered is plausible. Alternative dispersal routes, such as a corridor south of the Tibetan Plateau or directly across the plateau itself, are weakly supported or not supported at all by either reconstructions based on present-day predictors or projections onto paleoclimate.

This region may have acted as a dispersal corridor between East and Central Asia. The existence of such a corridor was likely not associated with the formation of continuous river basins, which are typically considered the main pathways for freshwater faunal dispersal and major biogeographic barriers (Starobogatov 1970; Graf 2007; Bolotov et al. 2020). However, freshwater mollusks, particularly members of the family Lymnaeidae, are capable of crossing such barriers by other means, for example, through passive transport by birds (Vinarski et al. 2017, 2021; Bespalaya et al. 2022).

The existence of such a “northern” corridor during the Pleistocene has been documented for certain terrestrial vertebrate taxa (Päckert et al. 2012; Choi et al. 2020), although no evidence supporting its existence for invertebrates or aquatic vertebrates was found. The most relevant comparison may be drawn with human dispersals, for which a dynamically opening and closing migration corridor between East and Central Asia has been proposed for the Middle and Late Pleistocene. For example, during Marine Isotope Stage 5 (ca. 130–70 ka), more humid and mild climatic conditions may have facilitated such dispersals (Li et al. 2019). The closure of this corridor during the Last Glacial Maximum (30–15 ka) is also consistent with the reduction in suitable habitat for Radix plicatula during this period.