The National Committee on Hydrology separates Thailand into 25 distinct hydrological units or river basins, which are used in this study as an established geographical reference. These units comprise the following rivers and drainage systems: Salween, Mekong, Kok, Shi, Moon, Ping, Wang, Yom, Nan, Chao Phraya, Sakaekrang, Pasak, Tha Chin, Mae Klong, Prachinburi, Bang Pakong, Tonle Sap, Peninsular East Coast, Phetchaburi, Peninsular West Coast, Southeast Coast, Tapi, Songkhla Lake, Pattani and Southwest Coast. These catchment and drainage systems are re-grouped here into seven areas, each with specific characteristics; refer to Figs 1, 4a, 8 for a graphical overview:

1 Central area: This is the most important area for Thailand, as it is an area without large water sources. The region, therefore, depends heavily on water from river basins upstream, such as Chao Phraya River as the main river of Thailand. The Chao Phraya begins at the confluence of the Ping and Nan rivers (Northern area) at Nakhon Sawan province. It flows from north to south from the central plains through Bangkok to the Gulf of Thailand.

2 Northern area: This area is a rich source of water for the central area (see above). For example, water of the Wang River flowing from north to south has its source in the Chiang Rai province. One of the principal cities along the river is Lampang, which is on the north bank of a curve in the river. From Lampang, the river flows southwards passing into Tak province. It joins the Ping River near Mae Salit north of the town of Tak. The Ping River originates in the Chiang Mai province, flowing through the provinces of Lamphun, Tak, and Kamphaeng Phet. The Nan River originates in the Nan province, subsequently draining the provinces Uttaradit, Phisanulok and Phichit. The Yom River joins the Nan River in the Chumsaeng district, Nakhon Sawam province. When the Nan River joins the Ping River it forms the Chao Phraya.

3 North-western area: This is a part of the drainage system of the Salween River, which flows into the neighbouring country of Myanmar.

4 Western area: This is part of the basin formed by the Me Klong River, which runs into the Gulf of Thailand.

5 North-eastern area: This is part of the Mekong river basin’s catchment area, which drains into the South China Sea.

6 Eastern area: An area characterized by many short rivers.

7 Southern area: Many short rivers and high annual rainfall characterize this area. There are a number of large water reservoirs.

Specimens of Tarebia granifera were collected throughout Thailand. For reference, we compare with samples available to MG from Timor Leste through the courtesy of Vince Kessner, who collected there recently. All samples were preserved in 95 % ethanol. Voucher specimens are kept in the collection of the Center of Natural History (CeNak), Zoological Museum, Universität Hamburg, Germany (ZMH) and the collection of the Parasitology and Medical Malacology Research Unit, Department of Biology, Faculty of Science, Silpakorn University, Thailand (SUT).

To reconstruct in detail the distributional range, in addition to own collecting activities in most parts of the region, material was analysed in several major museum collections, as well as literature records which were sufficiently verifiable as to the species identity (in general documented by descriptions and, even better, figures of shells collected).

Geographic coordinates of newly collected material were taken with a GPS device at the sampling site (WGS84 datum). Where GPS data for sampling sites were unavailable, coordinates were determined as accurately as possible from a map. Localities of the samples were mapped on a dot-by-dot basis on a public domain map (NaturalEarth, www.naturalearthdata.com) with ArcMap 10.4.1 (Esri Inc., Redlands, CA, USA). Final maps were compiled using Photoshop CS6 (Adobe Systems Inc., San José, CA, USA). The spelling of localities (whenever possible) follows GeoNames (http://www.geonames.org).

For climatic data, we used information from the climate of the world database (https://www.weatheronline.co.uk/reports/climate/Thailand).

The snails identified as belonging to Tarebia were grouped according to their morphological characteristics and geographic origin in four preliminary classes or morphs (see Fig. 2): 1.) specimens from Timor corresponding morphologically to T. granifera, i.e. without spiral pattern of narrow brown bands (Timor), 2.) specimens from Thailand corresponding to T. granifera (morph A), 3.) specimens from Thailand corresponding to T. lineata, i.e. specimens with a pattern of narrow brown spiral ridges (morph B), and 4.) Tarebia specimens from Thailand with broad brown spiral bands (morph C). In a second approach, specimens were grouped according to mitochondrial clades for morphological comparisons and analyses (see below).

The following biometrical parameters of the adult shells were taken with a digital calliper (accuracy: 0.1 mm): height of shell (h), width of shell (w), length of aperture (la), width of aperture (wa), height of body whorl (hbw), height of the last three whorls (l3w) and number of whorls (nw) (Fig. 3a). We were able to obtain these measurements from a total of 1,154 specimens. Analyses of shell parameters were performed using RStudio (RStudio Team 2016), with packages “ade4” (Chessel et al. 2004), “lawstat” (Hui et al. 2008), “agricolae” (Mendiburu 2010) “dunn.test” (Dinno 2017), and “car” (Fox and Weisberg 2011). For further testing the data set was partitioned in different predefined groups. These were: (a) four different morphs (see above) and (b) two different mitochondrial clades based on our molecular genetic analyses (see below). These different subdivisions of the data required slightly different approaches with regard to statistical testing. For the four different morphs, we first tested for normal distribution. Hence, we conducted the Shapiro-Wilk test for each subgroup individually. In case all measured variables for each group were normally distributed (Shapiro-Wilk-test p > 0.05), we performed an analysis of variance (ANOVA). If significant, it was followed by a Bonferroni-corrected LSD-Test. If at least one test for normal distribution (Shapiro-Wilk-test p < 0.05) was not significant, we instead deployed a Kruskal-Wallis-rank sum test. If the latter was found to be significant, a Bonferroni-corrected Dunn-test was conducted subsequently.

Figure 3. 

Biometrical parameters (a) and position of landmarks (b). Abbreviations: height of shell (h), width of shell (w), length of aperture (la), width of aperture (wa), height of body whorl (hbw) and height of last three whorls (l3w).

For the two different mitochondrial clades, the Shapiro-Wilk-test was performed on all measured variables for each group individually to test for normal distribution. If at least one group was not normally distributed, we conducted a Wilcoxon signed rank test with continuity correction, to test for significant differences between clades. If the data for both groups were normally distributed, a Levene-test based on absolute deviations from the mean was performed to check for homoscedasticity. In case homoscedasticity was detected, we tested different groups using a two-sample t-test. Otherwise, we performed Welch’s heteroscedastic t-test.

All available type specimens and the other examined material was photographed by remote shooting with EOS Utility 2.12.2.1 for Windows (Canon Inc., Tokyo, Japan) and Digital Photo Professional 3.12.51.2 for Windows (Canon Inc.) using a digital camera (EOS 5D MKII with Canon macro photograph lens MP-E 65 mm and compact macro lens EF 50 mm, Canon Inc.). Shell orientation was adjusted so that the apertural plane of the shell was perpendicular in relation to the optical axis of the camera and the shell’s columella parallel to the background. Photo stacks were assembled in Helicon Focus 5.3.14.2 for Windows (Helicon Soft Ltd., Kharkiv, Ukraine). The images were then edited with Photoshop CS6 (Adobe Systems Inc.).

A total of 1,169 standardized images of adult, unbroken shells could be included in our geometric morphometrics data set. Using tpsUtil version 1.74 (Rohlf 2017a), a tps-file including all specimens was assembled. We placed 15 landmarks (see Fig. 3b for landmark positions) with tpsDig2 version 2.30 (Rohlf 2017b). Data were analysed in RStudio, with “geomorph” (Adams and Otárola-Castillo 2013) and all packages listed for our biometry data analysis. After performing a Procrustes superimposition, a principal component analysis (PCA) was conducted to identify major axes of variance and to reduce dimensionality. Only axes with a relevant proportion of variance (> 0.05) were included. Following this procedure, we analysed the data set with partitions and tests as described above for the analyses biometric measurements (for each partition and principal component).

The content of the brood pouch was counted as best proxy for differences in the thiarid reproductive strategy following the method described in Glaubrecht et al. (2009) and Maaß and Glaubrecht (2012). The shells were cracked with a small vice, the operculum cut off from the posterior part of the foot using a scalpel and the soft body opened under a stereo microscope. After opening the brood pouch, which is located in the neck region of the female, care was taken to count all embryos and shelled juveniles contained within the marsupium, according to the nine standard size classes established for Thiaridae before by Glaubrecht et al. (2009): 1.) early embryos, 2.) late embryos, 3.) juveniles up to 0.5 mm, 4.) juveniles between 0.6 mm and 1.0 mm, 5.) juveniles between 1.1 mm and 1.5 mm, 6.) juveniles between 1.6 mm and 2.0 mm, 7.) juveniles between 2.1 mm and 2.5 mm, 8.) juveniles between 2.6 mm and 3.0 mm, and 9.) juveniles > 3.0 mm. We compared brood pouch contents for the different predefined morphological/geographic groups as described above, the main mitochondrial clades and for the different river systems in Thailand.

Sequences from 131 specimens of T. granifera from 95 populations in Thailand and 12 specimens from 11 populations in Timor Leste were generated (see Table 1). Two specimens of Thiara amarula (Linnaeus, 1748) were selected as outgroup. Total genomic DNA was extracted from ethanol-preserved foot tissue using a CTAB protocol (Winnepenninckx et al. 1993). For phylogenetic analyses fragments of the mitochondrial cytochrome c oxidase subunit 1 (cox1; 658 bp) gene using the primers LCO1490 (5’-GGT CAA CAA ATC ATA AAG ATA TTG G-3’, Folmer et al. 1994) and HCO2198var (5’-TAW ACT TCT GGG TGK CCA AAR AAT-3’, Rintelen et al. 2004) and the 16 S rRNA (16S; c. 780 bp aligned) gene using the primers 16S_F_Thia2 (5’-CTT YCG CAC TGA TGA TAG CTA G-3’, Rintelen, unpublished data, see Gimnich 2015) and H3059 (5’-CCG GTY TGA ACT CAG ATC ATG T-3’, Wilson et al. 2004) were amplified by PCR. Amplifications were conducted in 25 µl volumes containing, 2.5 µl 10× DreamTaq Green Buffer (Thermo Fisher Scientific, Waltham, MA, USA), 1.0 µl dNTP mix (5 mM each), 1.0 µl of each primer (10 µM), 0.2 µl of DreamTaq DNA polymerase (Thermo Fisher Scientific), 1.0 µl DNA template and 18.3 ddH₂O. After an initial denaturation step of 3 min at 94 °C, 35 cycles of 30 s at 94 °C, 60 s at 45–62 °C and 60–120 s at 72 °C were performed, followed by a final extension step of 5 min at 72 °C. Prior to sequencing, PCR products were enzymatically cleaned by adding 0.65 μl thermosensitive alkaline phosphatase (Thermo Fisher Scientific) and 0.35 μl exonuclease I (Thermo Fisher Scientific) to a 5 μl aliquot of the PCR reaction followed by an incubation step at 37°C for 15 min and enzyme inactivation at 85 °C for 15 min. Both strands of the amplified products were sequenced at Macrogen Europe Laboratory (Amsterdam, The Netherlands).

Forward and reverse strands were assembled using the program Geneious (Biomatters Limited, Auckland, New Zealand) and corrected by eye. The protein coding cox1 sequences were aligned with MUSCLE (Edgar 2004) as implemented in MEGA7 (Kumar et al. 2016) under default settings. The 16S sequences were aligned with MAFFT (Katoh and Standley 2013) using the Q-INS-i iterative refinement algorithm and otherwise default settings, because this algorithm has been described to perform better for the alignment of sequence data sets that may contain deletions and insertions than alternative multiple sequence alignment methods (Golubchik et al. 2007).

For information on vouchers and GenBank accession numbers, see Table 1. Pairwise genetic p-distances for the cox1 and 16S data sets were calculated with MEGA7.

Bayesian Inference (BI), Maximum likelihood (ML) and maximum parsimony (MP) approaches were used to reconstruct the phylogenetic relationships. The sequence data set was initially divided into four partitions for the nucleotide model-based ML and BI approaches: 1.) 1^st codon positions of cox1, 2.) 2^nd codon positions of cox1, 3.) 3^rd codon positions of cox1, and 4.) 16S. To select an appropriate partitioning scheme and/or evolutionary models for the mitochondrial sequences, the data set was analysed with PartitionFinder 2.1.1 (Lanfear et al. 2012) conducting an exhaustive search and allowing for separate estimation of branch lengths for each partition using the Bayesian information criterion as recommended by Luo et al. (2010) for model selection. Models to choose from were restricted to those available in MrBayes 3.2.6 (Ronquist et al. 2012) as well as in Garli 2.1 (Zwickl 2006). As best-fit partitioning scheme, the PartitionFinder analysis suggested to combine all predefined partitions into a single partition, with the HKY+G model as best-fit model under the Bayesian information criterion.

The BI analysis was performed using MrBayes 3.2.6. Metropolis-coupled Monte Carlo Markov chain (MC³) searches were run with four chains in two separate runs for 50,000,000 generations with default priors, trees and parameters sampled every 1,000 generations under default heating using the best-fit model as suggested by PartitionFinder. Diagnostic tools in MrBayes, including estimated sample size (ESS) values ≥ 200, were used to ensure that the MC³ searches had reached stationarity and convergence. The first 5,000,000 generations were discarded as burn-in.

Heuristic ML analysis was performed with Garli using the best-fit models as suggested by PartitionFinder. Support values were computed by bootstrapping with 1,000 replications.

Heuristic MP searches were carried out with PAUP v4.0b10 (Swofford 2002) using 100 random-addition-sequence replicates and TBR branch swapping. Support values were computed by bootstrapping with 1,000 replications.

Bayesian posterior probabilities (PP) values ≥ 0.95 and bootstrap (BS) values ≥ 70% and were interpreted as significant/meaningful support. BS values from the ML and MP analyses were mapped onto the Bayesian 50% majority-rule consensus tree with SumTrees 3.3.1, which is part of the Dendropy 3.8.0 package (Sukumaran and Holder 2010).

We used the General Mixed Yule-coalescent (Pons et al. 2006) in its Bayesian implementation (bGMYC) (Reid and Carstens 2012) and the Automatic Barcode Gap Discovery (ABGD) (Puillandre et al. 2012) with p-distances for DNA sequence-based species delimitation. The bGMYC method allows for taking phylogenetic uncertainty into account by basing the analysis on several ultrametric trees sampled from the same posterior distribution. We constructed ultrametric trees for the concatenated 16S and cox1 data set with Beast 2.4.1 (Bouckaert et al. 2014) assuming a strict clock and the same evolutionary model as in the Bayesian and ML analyses (root age was set to one using a lognormal prior). Chains were run for 10,000,000 generations discarding the first 50% of the generations as burn-in and sampling every 50,000^th tree resulting in a set of 100 ultrametric trees which were used in the bGMYC analyses. For each of the 100 ultrametric trees in the 16S and cox1 data set, the Markov-chain Monte Carlo sampler implemented in the bGMYC R package (Reid and Carstens 2012) was run for 100,000 generations, discarding the first 90,000 generations as burn-in and sampling every 100 generations.

We dated the divergence times for the main clades of Tarebia included in this study using the Bayesian algorithm implemented in Beast 2 based on the concatenated mitochondrial data assuming a strict molecular clock as the test implemented in MEGA 7 (Kumar et al. 2016) rejected a strict molecular clock at α = 0.05. The same partitioning scheme and nucleotide substitution models as in the MrBayes analysis described above were used. As tree prior the Yule speciation model was chosen. In the absences of fossil calibration points, a constant substitution rate of 1% per Ma was assumed as has previously been done by Köhler and Glaubrecht (2010) for related freshwater cerithioideans in the Pachychilidae. The Beast 2 analysis was run for 10,000,000 generations with a sampling frequency set to 10,000. Tracer 1.7 (Rambaut et al. 2018) was used to assess convergence of runs and to check whether effective sample sizes for all estimated parameters were above 200. A maximum clade credibility tree with median node heights was calculated with Treeannotator v1.8.2 from the BEAST 2 package discarding 10% of generations as burn-in.

In the present study, we examined phenotypically distinguishable shell morphs of yet another thiarid from Thailand, traditionally assigned to Tarebia granifera, in reference to samples from Timor Leste as known type locality of the nominal species, using biometry and geometric morphometrics in combination with phylogeographical analyses of molecular genetics and reproductive strategy.

We found Tarebia to be widespread in almost all freshwater bodies throughout Thailand, with a wide range of conchological variants or morphs, of which some closely resemble the types and topotypical material of granifera collected on Timor. While in Thailand Tarebia has been reported with only one species by Brandt (1974), distinct shell morphologies allow to distinguish phenotypically disparate morphs. Some of these have even been formally named as distinct species (albeit from other regions in Asia), based on ornamental features such as tubercles and/or nodules as well as the formation of elevated spiral ridges prominent in particular on the last body whorls. For example, the name lineata (as well as lateritia, traditionally used for Philippine forms) have frequently been applied to morphs and/or populations in the Oriental region. Subba Rao (1989) discussed that T. lineata was often synonymised with T. granifera, or treated as its variety (e.g. Benthem-Jutting 1959), although it is readily distinguished from the latter by the presence of the very distinct spiral ridges. Also, Appleton et al. (2009) for invading populations of T. granifera in South Africa described two distinct morphological variants found at different locations, among them also one with pronounced spiral ridges.

Applying a drainage-based phylogeographical as well as a biometrical approach, we were unable to find for the populations in Thailand any correlation of the morphs distinguished in this study based on discernable shell features as well as overall “Gestaltwahrnehmung” with any criteria deducible from our observations given above, neither with geographical occurrence or preferred habitat and substrate nor with the molecular genetic substructuring detected (see below). So, all available evidence points at the coexistence of different morphologies or disparate phenotypes in Tarebia granifera in this part of mainland Southeast Asia. However, in the absence of any of the discussed parameters or factors to be causally correlated with these morphological differences we are left with the hypothesis that they either qualify for reflecting phenotypical plasticity correlated with ecological variables in the habitat of the individual populations studied, and/or, alternatively, being correlated with the parthenogenetic reproduction discussed further below.

Biometric analyses are found useful tools for the study of characteristics that shape morphologically distinct entities, thus allowing to look into evolutionary pattern (e.g. Bocxlaer and Schultheiß 2010, Maaß and Glaubrecht 2012). Geometric analyses are used in addition to traditional morphometrics in order to compare in detail different populations and relationships among variable groupings (e.g. Rohlf and Marcus 1993, Sheets et al. 2006).

Although there are some differences in the biometric parameters and in the geometric morphometrics of Thai Tarebia, it is generally impossible to delimit distinct entities (in the sense of being at least indicative of the existence as biological species) based on these features, as all of them largely overlap (Figs 6, 7). Thus, our morphometric data do not support the distinction of T. lineata or other morphs from the nominal T. granifera, based on shell size and/or form. In addition, the same holds true for the two distinct molecular clades separated by mitochrondrial DNA sequences used in our study, for which we failed to find any diagnostic features in shell morphology or other phenotypical characteristics.

The measurement of shell height of T. granifera showed that they are within the size range previously reported as to vary between 6 to 44 mm (e.g. Abbott 1952, Brandt 1974). Also Bradstreet and Rogowsky (2012) reported on specimens of T. granifera to exhibit the same overall shape with an elongately or ovate-conoidal shell with the size index (L3W/W) in the order of 0.54–2.65 mm (Fig. 6d). Isnaningsih et al. (2017) found the shape of T. granifera from the Indonesian islands of Lombok, as well as from Banten and Maros, to be for the ratio of shell height to width 1.29–3.02 mm.

The results of geometric morphometrics revealed the overall shell shape of T. granifera from Thailand to be very similar to, and virtually undistinguishable from, conspecifics from Timor Leste (Fig. 6e). Thus, although T. granifera exhibit shell polymorphism this intra- and interpopulational variability in its shell characteristics does not allow for species-specific differentiation, as it was found, for example, in the thiarid Melanoides (e.g. Facon et al. 2003, Genner et al. 2004, Sorensen et al. 2005, Yousif et al. 2009).

In contrast to shell morphology (morphs A–C, or lineata vs. granifera phenotypes), we found based on molecular genetics strong indication as to the distinction of at least two natural entities within Tarebia in Thailand. As our analyses revealed, there is a most pronounced separation of two distinct mtDNA clades in this taxon, marked on the one hand by long branches in the resulting phylogenetic tree connecting these two clades, and on the other hand by very shorter branches within each of them (Figs 4, 5).

Therefore, our analyses would potentially allow for a more narrow species delimitiation within what has been to date traditionally treated in Thailand as T. granifera only (Brandt 1974). At the same time, the two clades correspond with a geographical separation into a northern and southern group. This is also reflected in ecology insofar, as both show a preference in altitude (see Fig. 4, insert). However, we propose that the latter reflects rather the occurrences in higher mountainous regions in the north than in the south of Thailand than a truly differential habitual preference. In contrast, the two genetically distinct lineages do neither match with features in shell morphology or biometry nor with differences in their reproductive strategies.

However, the p-distance of 13.8 % for cox1 and 10 % for 16S sequences has to be considered relatively high, hinting potentially at the existence of two genetically distinct species. However, a definite decision as to this species question in Tarebia in Thailand should remain open until the geographical distribution of genetically characterized populations of T. granifera and other congeneric forms is completely resolved and better understood within the entire autochthonous range in the Oriental region. Thus, it should be the privilege of a more comprehensive and in-depth analysis of the biogeographical situation based on an ongoing molecular genetic study (Glaubrecht unpubl. data).

While we found representatives of clade A in the northern tributaries of rivers such as the Chao Phraya and Mae Klong that run into the Gulf of Thailand, with only few others occurring at some localities in the south of Thailand (Figs 4, 8, 9a,b), those in clade B were found in the Salween River and the headwaters of Ping, Wang, Yom and Nan River. Accordingly, Tarebia snails from clade A are more frequent in the central to southern part of the country, whereas those from clade B are more frequent in the northern part. This overall geographic picture allows to attribute clade A as an element of the Sundaic region, given that it extends even further south and also comprises the Timor group (thus rendering it the nominal granifera), while clade B is mainly distributed in the Indochinese region (Figs 1, 4).

However, although being more frequent in the northern provinces, some representatives of clade B also occur in more southern locations, such as in the provinces Surat Thani (SUT 0516137), Nakhon Si Thammarat (SUT 0516139) and Phatthalung (SUT 0516138). We anticipate that this might reflect occurrences of passive dispersal, potentially via aquatic plant or other material or even transport by birds, rather than vicariance via the influence of sea level or tidal flows in drainage systems. The results of the median-joining haplotypes network and bGMYC analysis (Fig. 5) reveal that clade A and B exhibit many steps separating these two groups, and have low probabilities of conspecificity between clade A and B (p=0–0.05).

As we found in our molecular analyses this major split of clade A and B in Thai Tarebia to be as old as most likely c. 5.32 million years ago (Fig. 5d), it is worthwhile to look for a possible biogeographic explanation of the above distribution. In general, distinct faunal and floral assemblages are biogeographically restricted by barriers to dispersal such as characteristic geomorphological boundaries, even when individual taxa among each of the biota on either side often vary and may not all reflect the same discrete pattern. As Bruyn et al. (2005) pointed out, historical biogeography while providing crucial insights into the relationship between biological diversity and earth history, as a consequence has its limitations. However, patterns of intraspecific molecular variation may show unambiguous evidence for such historical divides, and can be used to test competing biogeographic hypotheses, such as e.g. the dispersal-vicariance debate, see e.g. Glaubrecht (2000) and Glaubrecht and Rintelen (2003) for limnic gastropods).

For the distributional pattern found in Tarebia in Thailand, a vicariant hypothesis can be formulated using a major biogeographic transition zone between the Sundaic and Indochinese biota, located just north of the Isthmus of Kra. It is interpreted as the result of Neogene marine transgressions that breached this isthmus in two locations for prolonged periods of time, i.e. more than 1 million year duration, as was shown e.g. by phylogeographic analyses of a freshwater decapod crustacean, the giant freshwater prawn Macrobrachium rosenbergii (cf. Bruyn et al. 2005). In his review Parnell (2013) examined, based on the relevant geological, geographical, climatic, biogeographic and sea-level data, the available evidence on the Isthmus of Kra as being a significant biogeographic divide on the Thai-Malay Peninsula and, thus, of mainland Southeast Asia. It is believed to be of the same scale as, e.g., the ‘Wallace’s Line’, albeit it remains less well-known and less well-studied, with its location and cause being still enigmatic. Dejtaradol et al. (2016) reported that population boundaries in birds did not coincide with the Isthmus of Kra, but instead were located north of the Thai‐Malay Peninsula in Central Thailand, while only one of four divides represented an Indochinese‐Sundaic transition. They supposed that different phylogeographical patterns among target species were presumably shaped by different ecological preferences in Pleistocene palaeohabitats. They found in bulbuls, as we suggest in analogy here for thiarid gastropods, Pliocene Indochinese‐Sundaic lineage divergence, for which they hypothesized that it coincides with strong vegetational changes on the Peninsula shaping two phytogeographical transitions. As distribution limits of bird species roughly coincide with these transition zones, the avifaunal Thai‐Malay transition represents apparently a broad zone rather than a sharp boundary.

While the separation of Tarebia and Thiara hint at a Late Miocene splitting event (anticipated to have occurred somewhere in the Indo-Malayan insular region of the Sunda and Sahul shelfs), our molecular and distributional data on T. granifera (Figs 4, 5d) suggest, with its two lineages in the north and south along the Thai peninsular mainland, to roughly correlate with a Late Miocene/Early Pliocene event (5.5–4.5 Mya). Thus, the separation of clade A and B can be hypothesized as resulting from a later marine transgression in the area to the north of today’s Isthmus of Kra that may have produced high sea-level stands with a seaway that dissected the Thai-Malay Peninsula for durations longer than one million years; see Bruyn et al (2005) and literature therein for further details as to the relevant geological data and discussion.

The fact that today the distributional boundaries of the two Tarebia populations in clade A and B do not coincide exactly with the position of the Isthmus of Kra, but are instead placed further to the north, could in this case be attributed to later palaeo-drainage differentiation in connection with orogenesis or other tectonic events in the mountainous central and northern regions of Thailand, as it was discussed using relevant geological and available biogeographical data, for example, from fishes and gastropods in Glaubrecht and Köhler (2004). Thus, although being today located north of the Thai‐Malay Peninsula in Central Thailand, the Isthmus of Kra and late Miocene/early Pliocene marine transgression might have caused in the freshwater thiarids of this region the separation of the Indochinese and Sundaic lineage within what has been regarded as Tarebia granifera to date.

Tarebia snails are all viviparous, i.e. they incubate embryos and later ontogenetic stages in an extra-uterine structure, called the subhemocoelic brood pouch, located at the back of the head in the female’s body running alongside and below, but being independent of the pallial organs (e.g. the gonoduct), and formed apparently by an invagination of the genital grove found in other oviparous cerithioidean gastropods (Glaubrecht 1996, 1999, 2006, 2011, Glaubrecht et al. 2009). Based on histology, for T. granifera Glaubrecht (1996) described an eu-viviparous strategy, involving matrotrophy (i.e. the nourishment by the female) of the progeny that develop in the subhemocoelic “marsupium” from early to late embryos and subsequently build their multi-whorled shells before hatching as crawling juveniles. This strategy, also known as typical for other thiarids such as e.g. Melanoides, is in contrast to an ovo-viviparous mode, reported e.g. for Thiara amarula and some other Australian thiarids, such as Stenomelania aspirans (see Schütt and Glaubrecht 1999, Glaubrecht et al. 2009, Maaß and Glaubrecht 2012).

In the Thai populations of Tarebia, as well as those from Timor, we found most if not all ontogenetic stages contained at the same time in the female’s marsupium, from early embryos to late embryos and shelled juveniles, in all morphs (A–C), both molecular genetic clades (A and B) and specimens from all drainage systems, without a clear-cut differentiation of this reproductive strategy. In particular, the ontogeny of T. granifera in Thailand is not obviously correlated to specific drainage systems, no matter where these water bodies eventually drain. Therefore, we conclude that Tarebia throughout its distributional range covered here is eu-viviparous, with only very few representatives in some populations (see Figs 8, 9, and more details above) that were found to only possess late and/or even early embryonic stages, respectively. This is in contrast to a pronounced correlation as to reproductive biology in Thiaridae from Australia (Glaubrecht et al. 2009), where all ovo-viviparous taxa that release veligers paradoxically exhibit very restricted distributional ranges in the Jardinian biogeographical region only. It is also in contrast to differences in the Thai thiarid Melanoides jugicostis (see Dechruksa et al. 2013), that was found to lack viviparous populations at least in some geographical regions and during some time of the year.

As in this later case, it could be hypothesized that any environmental factor might affect the reproductive strategy also in Tarebia. However, our analysis of representative climatic charts for the two parameters temperature and precipitation revealed no clear regional pattern of brood pouch content, as no correlation with the various ontogenetic stages were found across all locations in Thailand where T. granifera was sampled (see Figs 8, 9, 10, 11a). However, as reported above (see Figs 11b, 12a) some populations in rivers in the northwest (Pai, Moei, Ping), that were sampled essentially in the first half of the year (i.e. particularly early in the rainy season from April to June) exhibit a considerable amount of non-gravid specimens. The same might be true for some populations sampled during the early rainy season (April–July) in the Gulf of Thailand drainages, and to a lesser extent, too, in samples collected in May in the Andaman Sea drainages. It can be provisionally deduced from these data, that there might be a tendency for Tarebia females to be gravid especially during and after the end of the main or great rainy season in the second half of the year (and potentially during the dry season). In contrast to this temporal (spatial) hypothesis, we do not explain the frequency of non-gravid specimens as being indicative of the varying existence of males, as their occurrence would alternatively be regionally specific (in the northwest) and seasonal (little rainy season early in the year), which we doubt.

As in most (if not all) thiarids, Tarebia apparently lacks males in most populations, as we failed to find positive evidence for their existence. Parthenogenetic reproduction has gained much interest in the past in evolutionary biology, not only with respect to the origin of sex. Clonal reproduction in natural populations has obviously many advantages over sexual modes, with growth rates in the former often being much accelerated over the latter, as all individuals within the population are able to contribute (Maynard Smith 1978). In addition, these clones are considered instrumental in fast colonization of new habitats and areas, as even a single female can give rise to a new population (Baker 1955). Nevertheless, most faunas are dominated by sexually reproducing species, with asexual organisms being in the minority (Bell 1982).

Also in malacology there are some classical case studies, such as the New Zealand freshwater hydrobiid Potamopyrgus antipodarum (Jokela et al. 2003) or the thiarid Melanoides tuberculata (Jacob 1957, 1958, Berry and Kadri 1974, Ben-Ami and Heller 2005). However, in both cases reproduction is not exclusively parthenogenetic. In populations of Melanoides tuberculata, for example, the frequency of males was found to vary between 40 % in the French West Indies (Samadi et al. 1998) and up to 66 % in Israel (Livshits and Fishelson 1983, Heller and Farstey 1990).

It would be tempting to anticipate a similar phenomenon of T. granifera in Thailand and Timor Leste here from the varying frequencies (with up to 17.40 %) of non-gravid specimens. However, none ad hoc feature such as e.g. shell morphology between male and female could be differentiated in these aphallic cerithiodeans. So, in the present study we assumed not only any brood pouch-bearing snail to be female but also those without brood pouch as being non-gravid females rather than being rare males, for the reasons discussed above in connection with regional and/or climatic differences.

Given the prediction supported here that thiarid gastropods reproduce largely (if not completely) via parthenogenesis, the application in particular of the biological species concept is not made easy in case of thiarids. Morrison (1954) in discussing the enormous shell variability in thiarids in context with parthenogenesis, noted wisely that “wise indeed is the scientist who can tell whether a clone is a species or not, and be right every time, in the case of the Thiaridae”. This was shown, for example, for Melanoides tuberculata (Jacob 1957, 1958, Berry and Kadri 1974, Facon et al. 2003). Therefore, Stoddart (1985) preferred to apply instead of the biological species concept that of the evolutionary species (ESC) following Wiley (1978, 1981), as in his opinion this concept “stresses the relevance of the process of speciation to species definitions and provides the most appropriate framework for the taxonomy of asexual organisms”. Although this statement is debatable for several reasons, admittedly, the biological species concept (BSC) is also not without problems in application to Tarebia, as it explicitly uses the reproductive criterion in sexually reproducing organisms. The BSC was introduced by Mayr (1942) and since then widely discussed; see literature survey with references updated and discussion with respect to limnic gastropods e.g. in Glaubrecht (2004, 2009, 2010, 2011, Glaubrecht et al. 2009).

In case of the thiarids it remains to be seen in how far they are actually prone exclusively to parthenogenesis. For example, for populations of Melanoides tuberculata in Israel Ben-Ami and Heller (2005) reported sexually as well as asexually reproducing individuals, thus contradicting the general assumption that indeed all thiarids reproduce via apomixis. Apparently, at least in M. tuberculata there are both modes realized, securing the exchange of genetic information by sexual reproduction as was shown in earlier allozyme studies (Livshits and Fishelson 1983) and excluding the possibility of gynogenesis, i.e. parthenogenesis with the development of eggs to be induced by contact with sperm, though. Given the fact that we have (albeit indirect) evidence for the presence of males at least in low frequency, as is evident from published records, e.g. on Thiara amarula (see Healy and Glaubrecht 2018), as well as unpublished data, we here anticipate at least the occasional sexual reproduction in Thiaridae. As their species either maintain low levels of males in some populations, or by other means switch between asexual and sexual reproduction, there is no objection to not applying species concepts grounded on the reproductive criterion, as explicitly done under the BSC.