Morphological and DNA analyses reveal cryptic diversity in Anentome wykoffi (Brandt, 1974) (Gastropoda: Nassariidae), with descriptions of two new species from Thailand

Ethics

All animal experiments were undertaken according to guidelines approved by the Institute of Animals for Scientific Purpose Development (IAD), National Research Council of Thailand (permit number U1-07724-2561, issued to Nithinan Chomchoei). This study complied with all the relevant national regulations and institutional policies for the humane care and use of animals.

Collection abbreviations

CMUZ, Museum of Zoology, Chiang Mai University, Chiang Mai, Thailand.

MNHN, Muséum national ďHistoire naturelle, Paris, France.

SMF, Forschungsinstitut und Naturmuseum Senckenberg, Frankfurt am Main, Germany.

ZRC, Zoological Reference Collection, Lee Kong Chian Natural History Museum, National University of Singapore, Singapore.

Sampling and specimen identification

Anentome wykoffi specimens were collected in various river basins and drainage systems in Thailand (Fig. 1 and Table 1) using the count-per-minute method (Olivier and Schneiderman 1956). Thailand comprises five main river basins, including the Salween, Chao Phraya, Gulf of Thailand coast, Peninsular Malaysia and Mekong basins (Food and Agriculture Organization of the United Nations 2021) (Fig. 1), with the latter being further subdivided into three fish migration zones (upper, middle and lower Mekong; Fig. 1) (Mekong River Commission 2010). Specimens were preserved at −20°C until being used to study morphology and DNA.

Identification of the specimens was primarily based on a comparison with photographs of the type specimens (Holotype SMF 219824 and Paratype SMF 225241) and the original description (Brandt 1974). Anentome wykoffi was hitherto likely only illustrated by Brandt (1974). Yet the species can be distinguished from A. helena by the well-rounded whorls, narrower costulation, thinner texture, shorter siphonal canal and lower number of whorls (A. wykoffi = 5, A. helena = 8) (Brandt 1974). For A. wykoffi, characters such as carination of the shoulder of the whorls, number of spiral lines on the spira and pattern of the central cusps on the central tooth of the radula were used to distinguish putative morphotypes as follows:

Shell measurements and geometric morphometric analyses

Shell dimensions, including shell height (SH), shell width (SW), apertural height (AH), apertural width (AW), last whorl height (LWH), last whorl width (LWW) and siphonal canal width (SCW) (Fig. 2a) were measured in mm with a digital Vernier calliper (BEC, G14032021) following Lee et al. (2022). The shell height without the last whorl (SH-LWH) was calculated as shell height minus last whorl height and spira height (SPH) was calculated as shell height minus apertural height. The following ratios were calculated: SH/SW, AH/AW, SH/AH, SW/AW, SH/LWH, LWH/(SH–LWH), LWH/AH, SPH/AH, SW/LWW, LWW/AW, SW/SCW, AW/SCW, LWW/SCW, LWH/SH, AH/SH and SPH/SH. A one-way analysis of variance (ANOVA) was applied with IBM SPSS Statistics software (ver. 20, IBM Corp.) to analyse the variation in shell parameters between putative morphotypes.

Fig. 2.

Shell measurements (a) and landmarks (b) of A. wykoffi used in this study. Shell height (SH), shell width (SW), apertural height (AH), apertural width (AW), last whorl height (LWH), last whorl width (LWW), siphonal canal width (SCW), spira height (SPH) and shell height without the last whorl (SH–LWH)

Shells were photographed with a Nikon D5500 digital camera for the geometric morphometric study. The photos were randomly ordered in tpsUtil (ver. 1.82, see https://tpsutil.software.informer.com/; Rohlf 2015). Fifteen landmarks on the shells were digitised on photos by the same person using the program tpsDig2 (ver. 2.29, see https://tpsdig2.software.informer.com/; Rohlf 2015) (Fig. 2b). Geometric morphometric analyses were performed with the software package MORPHOJ (ver. 1.06d, see https://morphometrics.uk/MorphoJ_guide/frameset.htm?changelog.htm; Klingenberg 2011). Canonical variance analysis (CVA) based on shell shape variation between putative morphotypes was examined. Pairwise Mahalanobis distances and Procrustes distances between putative morphotypes were assessed for significant differences by a permutation test (10 000 iterations). Photos of the type specimens from this study were deposited in Morphobank (http://morphobank.org/permalink/?P4892).

Scanning electron microscopic (SEM) study of the shell and radula

Shells were washed with distilled water several times and air dried to study shell microsculpture (Coppois and De Vos 1986). The shells subsequently were mounted on aluminium stubs using carbon conductive adhesive tape. The mounted shells were coated with gold and studied with an LV-Scanning Electron Microscope (JSM 5910 LV) at the Electron Microscope Research and Service Center, Faculty of Science, Chiang Mai University.

Radulae were extracted under a stereomicroscope, boiled with 10% sodium hydroxide for 5 min, washed with distilled water and observed under a stereomicroscope. Radulae were dehydrated by immersion in increasing concentrations of ethanol (10, 20, 30, 50, 70, 85 and 95%) before SEM, with each step taking 5 min (Yang and Zhang 2011). Radulae were mounted on aluminium stubs using carbon conductive adhesive tape, coated with gold and observed with an LV-Scanning Electron Microscope (JSM 5910 LV) as listed above.

DNA extraction, amplification and sequencing

Genomic DNA of 21 specimens of Anentome aff. wykoffi was extracted from a small piece of the foot muscle tissue using 150 µL of 5% of Chelex 100 and 3 µL of Proteinase K. Samples were incubated at 55°C for 1 h, 95°C for 30 min and centrifuged at 11 336g at room temperature for 1 min. DNA extracted from the supernatant was stored at −20°C until use. A fragment of the mitochondrial cytochrome c oxidase subunit I (COI) gene and a fragment of the nuclear 28S rRNA gene were amplified by PCR using the primers LCOI490 (5′-GGT CAA CAA ATC ATA AAG ATA TTG G-3′) and HCO2198 (5′-TAA ACT TCA GGG TGA CCA AAA AAT CA-3′) for COI (Folmer et al. 1994), and C1 (5′-ACC CGC TGA ATT TAA GCA T-3′) and D2 (5′-TCC GTG TTT CAA GAC GG-3′) for 28S rRNA (Chisholm et al. 2001). PCR amplification was performed in reaction volumes of 20 µL containing 10 µL of 2× Illustra Hot Start Master Mix (GE Healthcare), 1 µL of each primer (10 mM), ~1 ng of DNA template and deionised water up to 20 µL. For COI, thermal cycling involved 94°C for 2 min, followed by 36 cycles of 94°C for 30 s, 42°C for 2 min and 72°C for 2 min, and a final extension step of 72°C for 5 min. For 28S rRNA, thermal cycling involved one cycle of 95°C for 3 min, 45°C for 2 min and 72°C for 90 s, followed by four cycles of 95°C for 45 s, 50°C for 45 s and 72°C for 90 s, in turn followed by 25 cycles of 95°C for 20 s, 52°C for 20 s and 72°C for 90 s and a final extension step of 72°C for 5 min. The amplicon size was checked by 1% (w/v) agarose gel electrophoresis using 1× TBE buffer. Gels were run at 100 V for 20 min and visualised with RedSafe nucleic acid staining solution and UV transillumination. The PCR products were purified using the BigDye Terminator (ver. 3.1) cycle sequencing kit chemistry. The amplified PCR products were directly cycle-sequenced using the original amplification primers with the sequencing reaction products run on 1st BASE DNA Sequencing Services (Applied Biosystems).

Phylogenetic analyses of DNA sequences

DNA sequences were edited and aligned with ClustalW (ver. 2, see http://www.clustal.org/clustal2/; Thompson et al. 1994) as implemented in MEGA (ver. 7, see https://www.megasoftware.net/; Kumar et al. 2016) and improved manually. COI sequences were checked for stop codons and frameshift mutations. 28S rRNA sequences were used after complete deletion of missing-information and alignment gaps (indels). All new sequences were deposited in GenBank (Table 1).

Phylogenetic trees were inferred using neighbour-joining (NJ) (Saitou and Nei 1987), maximum likelihood (ML) (Sullivan 2005) and Bayesian inference (BI) (Ronquist et al. 2012). The genus Nassarius was used as an outgroup. The best-fit evolutionary substitution model based on the Akaike Information Criterion (Akaike 1974) as implemented in jModeltest2 (ver. 2.1.10, see https://github.com/ddarriba/jmodeltest2; Darriba et al. 2012) was GTR + I + G for the tree datasets: COI, 28S rRNA and the concatenated COI–28S alignment. NJ- trees were constructed and pairwise genetic distances based on the COI sequences were calculated using the Kimura 2-parameter model (K2P) (Kimura 1980) in MEGA (ver. 7; Kumar et al. 2016) with 1000 bootstrap replicates. ML-trees were constructed using PhyML (ver. 3, see http://www.atgc-montpellier.fr/phyml/; Guindon et al. 2009) with 1000 bootstrap replicates. Bootstrap values higher than 70% were considered to provide strong support (Hillis and Bull 1993). Bayesian inference was performed using MrBayes (ver. 3.2.7, see https://nbisweden.github.io/MrBayes/index.html; Ronquist et al. 2012). The Markov Chain Monte Carlo (MCMC) search was run with four chains for 7 000 000 generations with the heating parameter set at 0.09, tree sampling every 100 generations and burn-in set at 25%. Posterior probabilities were considered significant when ≥0.95 (San Mauro and Agorreta 2010). Tree topologies were drawn with FigTree (ver. 1.4.3, A. Rambaut, see http://tree.bio.ed.ac.uk/software/figtree/).

Haplotype network analyses

COI haplotypes were analysed using DnaSP (ver. 6.12.03, see http://www.ub.edu/dnasp/; Rozas et al. 2017). A haplotype network was generated and visualised with the TCS method (Clement et al. 2000) as implemented in PopART (ver. 1.7, see https://popart.maths.otago.ac.nz/; Leigh and Bryant 2015). River drainage systems where the detected haplotypes of A. wykoffi were found (Table 1) were used as traits and employed in building the haplotype network (Yodsiri et al. 2017).

Species delimitation

Species delimitation methods were applied to the COI and 28S rRNA data using the (i) Assemble Species by Automatic Partitioning (ASAP) (Puillandre et al. 2021), (ii) generalised mixed Yule-Coalescent (GMYC) (Fujisawa and Barraclough 2013) and (iii) Bayesian Poisson tree processes (bPTP) (Zhang et al. 2013) methods.

ASAP was run on the web server (see https://bioinfo.mnhn.fr/abi/public/asap/) with the default settings and the K2P (Kimura 1980) substitution model.

The GMYC and bPTP methods require the input of an ultrametric tree from the BEAST package (ver. 1.10.4, see https://beast.community/; Suchard et al. 2018). An XML file was made with BEAUti (ver. 1.10.4, see https://beast.community/beauti) using the relaxed log-normal clock algorithm with the GTR + I + G model. The tree was reconstructed for 1 × 10⁷ generations with sampling every 1000 steps. The output tree was analysed with TreeAnnotator (ver. 1.10.4, see https://beast.community/treeannotator). The tree file was displayed in FigTree (ver. 1.4.3). GMYC was performed using both a single and multiple thresholds, and run on the web server (see https://species.h-its.org/gmyc/). bPTP was carried out on the bPTP web server (see https://species.h-its.org/ptp/).

Sampling

In total, 125 specimens of Anentome aff. wykoffi were collected from 12 localities (Table 1, Fig. 1) in Thailand. The species seems to occur in at least five drainage systems in three river basins, viz. the Wang drainage in the Chao Phraya basin, the Chi, Mun and Mekong drainages in the middle Mekong basin, and the Tonle Sap drainage in the lower Mekong basin (Fig. 1). Based on the intraspecific variation we observed in A. wykoffi, we defined three putative morphotypes that could be distinguished by combining three characters: the presence or absence of a carinated shoulder, the number of spiral lines on the spira and the pattern of the central cusps on the central tooth of the radula.

Morphological analysis

The one-way ANOVA (with Tukey’s post hoc test for multiple comparisons) revealed that the shells of the three putative morphotypes differed significantly (P < 0.01) from each other for 19 out of the 25 shell parameters studied (not AH, SH/SW, SW/AW, LWW/AW, SW/SCW and AW/SCW; Table 2).

In the geometric morphometric analysis, the CVA based on shell shape variation with the three putative morphotypes of A. wykoffi as a priori groups provided a graphic display of the shape differences based on the two CV variables obtained (Fig. 3a). CV1 explained 76.90% of the total variance in shape variation, with an eigenvalue of 0.9684 and CV2 explained the remaining 23.10% of the total variance in shape variation, with an eigenvalue of 0.2909. Overall, the shell shapes of the three putative morphotypes were rather well distinguished, even if these showed some overlap. This distinction was confirmed by the significant differences among the putative morphotypes as shown by the permutation tests of the Mahalanobis and Procrustes distances (P < 0.0001).

Fig. 3.

Geometric morphometric study of shell shape variation in A. wykoffi. (a) Plots of individual scores for the canonical variates (CVs) that were derived from the canonical variate analysis (CVA); blue, green and red circles represent individuals with morphotypes a, b and c respectively, and black squares represent the type specimens of A. wykoffi (Holotype SMF 219824 and Paratype SMF 225241). (b) Wireframe modifications derived from the CVA showing the shell shape changes (black solid line) from the consensus configuration (red dashed line) in relation to each of the negative and positive portions of the CVs. Shape changes along CV1 are shown on the left and those along CV2 on the right.

The wireframe modifications in Fig. 3b describe the shell shape changes from the consensus configuration along the two CVs. Individuals located in the positive portion of the axis for CV1 had a narrower shell than individuals located in the negative portion. Individuals located in the negative portion of the axis for CV2 had a wider aperture than individuals located in the positive portion.

Phylogenetic analyses of DNA sequences

Nucleotide sequences of the COI (657 bp), 28S rRNA (728 bp) and concatenated dataset of the two fragments (1385 bp) were aligned, along with Nassarius arcularia and N. boucheti as outgroup taxa. Not unexpectedly, the trees based on the separate analyses of COI and 28S rRNA were somewhat less well resolved but remained largely congruent or did not significantly contradict the concatenated analyses.

The NJ, ML and BI trees based on the concatenated dataset revealed that the four species of the genus Anentome (A. cambojiensis, A. helena, A. costulata and A. wykoffi) form a maximally supported clade in all analyses (Fig. 4).

Fig. 4.

BI tree of A. wykoffi and related species based on 1385 nucleotide sites in a concatenated COI (657 bp) and 28S rRNA (728 bp) alignment. Numbers at the branches are NJ bootstrap (BP)/ML bootstrap (BP)/BI posterior probability (PP) values. The map shows the sampling localities (numbered from 1 to 12 in both the map and the phylogenetic tree) with their morphotypes (clade colours).

The three morphotypes of A. wykoffi formed three well-supported clades, making A. wykoffi paraphyletic with the inclusion of A. cambojiensis (Fig. 4). Clade A (morphotype A) included mainly specimens from the Chi, Mun and Mekong drainages (localities 2–9). Clade B (morphotype B) included specimens from the Tonle Sap (locality 12) and Mun drainages (locality 10–11). Clade C (morphotype C) included specimens from the Wang drainage (locality 1) and formed a well-supported clade with A. cambojiensis.

The COI K2P-distances among Anentome species ranged from 3.65 ± 0.82 to 15.82 ± 1.79%. The K2P-distances among three clades of A. wykoffi ranged from 3.74 ± 0.75 to 7.17 ± 0.98%. The intraspecific K2P-distances within A. wykoffi clades A, B and C were 1.82 ± 0.31, 0.10 ± 0.13 and 0.00 ± 0.00% respectively (Table 3).

Haplotype network analysis

The TCS haplotype network of COI sequences (Fig. 5) suggests three haplogroups that correspond with the three putative morphotypes A, B and C. Haplogroup A comprises 11 haplotypes (Hap1-11) from the Chi, Mun and Mekong drainages (localities 2–9), with hap1 and hap2 being somewhat differentiated from the other A group haplotypes by 11 mutational steps. Haplogroup A differs from haplogroups B and C by 27 and 41 mutational steps respectively. Haplogroup B comprises hap12 (Mun drainage; localities 10 and 11) and hap13 (Tonle Sap drainage; locality 12) that are separated by only two mutational steps. Haplogroup C only comprises hap14 (Wang drainage; locality 1) that differs from haplogroup B in 28 mutational steps.

Fig. 5.

(a) Distribution of the haplotypes of A. wykoffi in the different drainage systems in Thailand (the numbers in the map are the locality numbers). (b) TCS network of COI haplotypes of A. wykoffi. Drainage systems are indicated by different colours. Circles of different sizes represent the number of individuals with a given haplotype. Hatch marks along the branches indicate the numbers of mutational steps and black dots represent hypothetical haplotypes (Hap, haplotype).

Species delimitation

The three species delimitation methods (ASAP, GMYC and bPTP) based on COI sequences consistently suggested that the three clades originally identified as A. wykoffi represent putative species, coinciding with morphotypes A, B and C. These are supported by morphological characteristics and can be distinguished by the carinated shoulder, number of spiral lines on the spira and pattern of the central cusps on the central tooth of the radula (Fig. 6). However, both GMYC and bPTP further subdivided clade A into three putative species (A-1, A-2 and A-3; Fig. 6). Yet these three putative species showed very limited COI sequence divergence (1.82%) and appeared morphologically indistinguishable, therefore we did not retain these for taxonomic treatment. Besides, ASAP and GMYC for 28S rRNA sequences failed to separate the three clades of A. wykoffi into three putative species. Conversely, the three morphotypes and COI–28S rRNA clades are proposed to be different species, two of which are new to science.

Fig. 6.

Species delimitation analyses of A. wykoffi based on COI and 28S rRNA sequences using ASAP, GMYC and bPTP. The putative species are labelled with letters and numbers. Posterior probabilities of the BI tree are provided at the branches. The upper left part of the figure provides a synopsis of three morphological characters by which the three morphotypes are primarily defined: (1) carinated shoulder, (2) number of spiral lines on spira, and (3) pattern of the central cusps on the central tooth of the radula.

Family NASSARIIDAE Iredale, 1916 (1835)

Genus Anentome Cossmann, 1901

Anentome wykoffi (Brandt, 1974) (sensu stricto)

(Fig. 4, clade A, and 7a, b, e, h.)

ZooBank: urn:lsid:zoobank.org:act:3307039C-6327-42CB-9E62-A2E600873001

Clea (Anentome) wykoffi Brandt, 1974, p. 204, pl. 15, fig. 69. Type locality: Mekong at Bandan in Thailand.

Clea wykoffi (Brandt, 1974) – Kittivorachate and Yangyuen (2004), p. 133, Ubol Ratana Dam, Khon Kaen, Thailand.

Anentome wykoffi (Brandt, 1974) – Chomchoei et al. (2022), table 1, fig. 1, Ban Kong, Nong Ruea, Khon Kaen, Thailand; Khemmarat, Ubon Ratchathani, Thailand; Khong Chiam, Ubon Ratchathani, Thailand.

Material examined

Type material

Holotype. SMF 219824 (Fig. 7a) and Paratype SMF 225241 from Mekong at Bandan in Thailand. Non-type material: THAILAND: Ubol Ratana Dam, Khon Kaen: CMUZ A213003 (1 shell) and CMUZ A213031 (3 specimens in −20°C), Huai Kut Daeng Reservoir, Roi Ed: CMUZ A193004 (1 shell), CMUZ A193005 (1 shell) and CMUZ A193032 (5 specimens in −20°C), Phaya Thaen Public Park, Yasothorn: CMUZ A193006 (1 shell), CMUZ A193007 (1 shell), CMUZ A193029 (1 shell), CMUZ A193030 (1 shell) and CMUZ A193033 (8 specimens in −20°C), Mekong River, Khemmarat, Ubon Ratchathani: CMUZ A213013 (1 shell) and CMUZ A213034 (8 specimens in −20°C), Wat Ban Dan, Ubon Ratchathani: CMUZ A213014 (1 shell), CMUZ A213015 (1 shell; Fig. 7b), CMUZ A213016 (1 shell) and CMUZ A213035 (7 specimens in −20°C), Mun River Bridge, Khong Chiam, Ubon Ratchathani: CMUZ A213008 (1 shell), CMUZ A213009 (1 shell) and CMUZ A213036 (8 specimens in −20°C), Sirindhorn Dam, Ubon Ratchathani: CMUZ A193010 (1 shell), CMUZ A193011 (1 shell) and CMUZ A193037 (11 specimens in −20°C), Lam Dom Noi, Ubon Ratchathani: CMUZ A213012 (1 shell) and CMUZ A213038 (2 specimens in −20°C).

Fig. 7.

Morphological comparison of species formerly grouped under the name Anentome wykoffi. (a) A. wykoffi holotype SMF 219824, (b) A. wykoffi CMUZ A213015 (Clade A in Fig. 4), (c) A. longispira sp. nov. (Clade B in Fig. 4) holotype CMUZ A213021, and (d) A. khelangensis sp. nov. (Clade C in Fig. 4) holotype CMUZ A193001. SEM images of the shell apex (e–g) and radula (h–j) of A. wykoffi (clade A) (e, h), A. longispira sp. nov. (clade B) (f, i) and A. khelangensis sp. nov. (clade C) (g, j) (white arrow: the carinated shoulder present; red arrow: the carinated shoulder absent).

Anentome longispira Chomchoei & Nantarat, sp. nov.

(Fig. 4, clade B, and 7c, f, i.)

ZooBank: urn:lsid:zoobank.org:act:7FF63177-8938-4E8C-8085-08A4D8CA6BE8

Anentome wykoffi (Brandt, 1974) – Chomchoei et al. (2022), table 1, fig. 1, Dan, Kap Choeng, Surin, Thailand; Non Din Daeng, Buri Ram, Thailand; Aranyaprathet, Aranyaprathet, Sa Kaeo, Thailand.

Anentome khelangensis Chomchoei & Nantarat, sp. nov.

(Fig. 4, clade C, and 7d, g, j.)

ZooBank: urn:lsid:zoobank.org:act:72E73AF2-5CEA-44DD-A6CA-E5A90F798DD0