- Open Access
Re-evaluating the case for poecilogony in the gastropod Planaxis sulcatus (Cerithioidea, Planaxidae)
BMC Ecology and Evolution volume 22, Article number: 13 (2022)
Planaxis sulcatus has been touted as a textbook example of poecilogony, with members of this wide-ranging Indo-Pacific marine gastropod said to produce free-swimming veligers as well as brooded juveniles. A recent paper by Wiggering et al. (BMC Evol Biol 20:76, 2020) assessed a mitochondrial gene phylogeny based on partial COI and 16S rRNA sequences for 31 individuals supplemented by observations from the brood pouch of 64 mostly unsequenced individuals. ABGD and bGYMC supported three reciprocally monophyletic clades, with two distributed in the Indo-Pacific, and one restricted to the northern Indian Ocean and Red Sea. Given an apparent lack of correlation between clade membership and morphological differentiation or mode of development, the reported 3.08% maximum K2P model-corrected genetic divergence in COI among all specimens was concluded to represent population structuring. Hence, the hypothesis that phylogenetic structure is evidence of cryptic species was rejected and P. sulcatus was concluded to represent a case of geographic poecilogony.
Our goal was to reassess the case for poecilogony in Planaxis sulcatus with a larger molecular dataset and expanded geographic coverage. We sequenced an additional 55 individuals and included published and unpublished sequence data from other sources, including from Wiggering et al. Our dataset comprised 108 individuals (88 COI, 81 16S rRNA) and included nine countries unrepresented in the previous study. The expanded molecular dataset yielded a maximum K2P model-corrected genetic divergence among all sequenced specimens of 12.09%. The value of 3.08% erroneously reported by Wiggering et al. is the prior maximal distance value that yields a single-species partition in ABGD, and not the maximum K2P intraspecific divergence that can be calculated for the dataset. The bGMYC analysis recognized between two and six subdivisions, while the best-scoring ASAP partitions recognized two, four, or five subdivisions, not all of which were robustly supported in Bayesian and maximum likelihood phylogenetic analyses of the concatenated and single gene datasets. These hypotheses yielded maximum intra-clade genetic distances in COI of 2.56–6.19%, which are more consistent with hypothesized species-level thresholds for marine caenogastropods.
Based on our analyses of a more comprehensive dataset, we conclude that the evidence marshalled by Wiggering et al. in support of Planaxis sulcatus comprising a single widespread, highly variable species with geographic poecilogony is unconvincing and requires further investigation in an integrative taxonomic framework.
Marine invertebrates can undergo one of two types of larval development . In the planktotrophic type of development, a free-swimming planktonic larva feeds on phytoplankton, typically during a period of several weeks or months, or occasionally years; in the non-planktotrophic (also called lecithotrophic or, sometimes incorrectly, direct) type of development, the larva reaches metamorphosis without phytoplanktonic food uptake  and typically spends a few days or less in the water column. The mode of larval development is a species attribute that is usually difficult to determine as, in most phyla, it requires collecting eggs from spawning adults and raising the larvae. In shelled Gastropoda, however, as a consequence of the accretionary growth of the shell, the protoconch—or larval shell, very frequently retained at the apex of the adult shell—reflects the mode of development; a so-called multispiral protoconch is indicative of planktotrophic larval development, and a so-called paucispiral protoconch is indicative of non-planktotrophic larval development . This extraordinary correlation can be traced in the fossil record , and gastropod protoconchs can be attributed to either a planktotrophic or a non-planktotrophic mode of development already in Palaeozoic fossils. Not surprisingly, the protoconch has been used as a species-specific character, with many examples of species pairs diverging in inferred development type, and developmental shifts suggested to represent a major driver of speciation (see e.g. ).
The literature contains reports of rare cases when larvae with the two modes of development are produced by the same individual, population, or species [6, 7], a phenomenon termed poecilogony. Notably, even among the rare and established cases of poecilogony, there is still only a single species (Alderia willowi Krug, Ellingson, Burton & Valdés, 2007) in which one individual has been shown to vary the development mode of its offspring . However, if the phenomenon were widespread, this would of course severely weaken the value of the multispiral/paucispiral protoconch dichotomy as a taxonomic character. Whereas poecilogony has been confirmed in several species of sea-slugs [9,10,11]—and, outside molluscs, in polychaetes —, alleged cases of poecilogony in shelled gastropods were reviewed in pre-molecular times by Hoagland and Robertson  and Bouchet  who independently concluded that there were no definitive cases of poecilogony in marine shelled gastropods. However, McDonald et al.  and Russini et al.  did provide molecular data suggesting the existence of very rare cases of poecilogony in Caenogastropoda, although these represented different taxa from the cases traditionally reported in the historical literature.
One such historically reported case concerns the marine snail Planaxis sulcatus (Born, 1778), a common intertidal gastropod with a broad Indo-Pacific distribution, from the East African coasts of the Indian Ocean and the satellite Red Sea and Persian Gulf, to the western central Pacific. Risbec  and Thorson  reported contradictory observations on its mode of development in New Caledonia and the Persian Gulf, respectively. In New Caledonia, Risbec described a classical planktotrophic larval development, whereas in the Persian Gulf the same species was reported to incubate the young, which feed on nurse eggs and metamorphose before hatching. Although Bouchet  hypothesized the existence of two cryptic species, each with its own species-specific mode of development, the Planaxis sulcatus case has been highlighted by Wiggering et al.  as a “textbook example for poecilogony”.
Planaxidae constitute a small family of tropical/temperate cerithioidean gastropods including about 60 Recent species. They are all restricted to the upper/middle intertidal, where they may form large aggregations, either on/under stones (Planaxinae) or in crevices (Fossarinae), where they presumably graze on the biofilm. Many species brood their fertilized eggs in a subhaemocoelic brood pouch located in the neck region of the headfoot, until larvae or juveniles are released.
Wiggering et al.  observed that different populations of Planaxis sulcatus exhibit different modes of larval development, with Western Indian Ocean and Red Sea populations releasing large, shelled juveniles from the brood pouch, whereas the Indo-West Pacific populations release planktotrophic veliger larvae. They sequenced 31 specimens of P. sulcatus from throughout the geographical range and concluded that their data confirmed P. sulcatus as a single species rather than a group of cryptic species, thus representing a case of geographic poecilogony.
If poecilogony occurs in Planaxis sulcatus, and other caenogastropods, then the paradigm of protoconch dichotomy (multispiral v. paucispiral) as a species-specific character would be severely weakened, with profound consequences for the alpha-taxonomy of extant and fossil molluscs. The genetic mechanisms involved in larval ecology of marine invertebrates (including intraspecific variation) are still largely unknown, and poecilogonous taxa would be the best models for such studies . However, the acceptance of this extraordinary claim requires that the evidence and its interpretation are unambiguous: for poecilogony to be robustly inferred, the two modes of larval development must be found to occur—or at least be reliably inferred to occur—among the offspring of the same individual, or among individuals of the same or different populations of a single species.
The objective of the present work is thus to review the case for alleged poecilogony in Planaxis sulcatus with an expanded molecular dataset, and to test the alternative hypothesis that P. sulcatus in fact comprises a complex of cryptic species potentially differing in development.
Wiggering et al.  generated a dataset of partial COI and 16S rRNA mitochondrial gene sequences from 31 museum specimens. Owing to difficulties in extracting and amplifying DNA from historical material, they were successful in obtaining only 16 COI and 28 16S rRNA sequences. We were able to augment this from published and unpublished sources, and with newly generated sequences from recently collected specimens in the collections of the Muséum national d’Histoire naturelle in Paris (MNHN), National Museum of Natural History in Washington, DC (USNM), and Florida Museum of Natural History (UF). Our dataset comprises 108 specimens of Planaxis sulcatus from three biogeographic realms (Western Indo-Pacific, Central Indo-Pacific, and Temperate Northern Pacific), plus one specimen collected in the Eastern Indo-Pacific (Palmyra Atoll) (following ; Fig. 1). Sequences for 53 specimens were retrieved from GenBank and BOLD, including those published by Wiggering et al. , and unpublished sequences for seven individuals were kindly shared with us by colleagues. Sequences for the two outgroups used by Wiggering et al. , Supplanaxis niger and Planaxis planicostatus, and for Planaxis sp. were also downloaded from GenBank. See Table 1 for details.
Molecular analyses and sequence alignment
Thirty-nine of the 55 specimens sequenced for this study were obtained during expeditions organized by the MNHN and Pro-Natura International as part of the Our Planet Reviewed program. These specimens were anesthetized using magnesium chloride (MgCl2) or were microwaved to separate the animal from the shell . Tissue clips of foot tissue were preserved in 95–98% ethanol. Specimens in the USNM were cracked and preserved whole in 95% ethanol.
Two labs (USNM; Sapienza University of Rome) contributed sequences for this study using slightly different protocols. At the USNM, whole genomic DNA was extracted from a ~ 1 mm3 tissue clip of the foot using an Autogenprep965 (Autogen, Holliston, MA) automated phenol:chloroform extraction with a final elution volume of 50 µL. A 691 base pair (bp) fragment of cytochrome c oxidase subunit I (COI) was amplified using the jgLCOI primer  in combination with Cerithioid_COIR ; a 505–508 bp fragment of 16S ribosomal DNA was amplified with the universal 16S AR/BR primers . PCR reactions were performed with 1 µL of undiluted DNA template in 20 µL reactions. Reaction volumes for COI consisted of 10 µL of Promega Go-Taq Hotstart Master Mix, 0.15 µM each primer, 0.25 µg/µL BSA, 1.25% DMSO and an amplification regime of an initial denaturation at 95 °C for 7 min, followed by 45 cycles of denaturation at 95 °C for 45 s, annealing at 42 °C for 45 s, extension at 72 °C for 1 min and a final extension at 72 °C for 3 min. Reaction volumes for 16S were 1x Biolase (Bioline, Taunton, MA) reaction buffer, 500 µM dNTPs, 3 mM MgCl2, 0.15 µM each primer, 0.25 µg/µL BSA, 1 unit Biolase DNA polymerase and an amplification regime of initial denaturation at 95 °C for 5 min, followed by 35 cycles of denaturation at 95 °C for 30 s, annealing at 48 °C for 30 s and extension at 72 °C for 45 s, followed by a final extension at 72 °C for 5 min. PCR products were purified using the ExoSAP-IT protocol (USB Corporation). BigDye 3.1 (ABI, Foster City, CA) sequencing reactions and sequencing on an ABI 3730XL DNA analyzer capillary array were done following manufacturer’s instructions. At Sapienza University of Rome, whole genomic DNA was extracted from a ~ 1 mm3 tissue clip of foot tissue by using a ‘salting out’ protocol , with a final elution of 50 µL. A 658 bp of COI was amplified using the jgLCOI and jgHCO primers ; a ~ 800 bp fragment of 16S ribosomal DNA was amplified with the 16SA  and CGLeuR  primers. PCR reactions were performed with 1 µL of undiluted DNA template in 25 µL reactions. Reaction volumes consisted of 2.5 µL of 10x NH4 Reaction Buffer, 2.5 µL of 50 mM MgCl2 Solution, 0.15 µL of BIOTAQ DNA Polymerase, 0.4 µL of each 25 pM primer solution, 1 µL of 10% BSA solution, 0.5 µL of 10 mM nucleotide mix solution. PCR conditions for COI followed a “touchdown” profile as in , while for 16S were as follows: initial denaturation (94 °C/4′); 35 cycles of denaturation (94 °C/30''), annealing (52 °C/40″), and extension (94 °C/1′); final extension (72 °C/10′). PCR products were purified using ExoSAP-IT (USB Corporation) and sequenced at Macrogen, Inc.
Amplicons were sequenced in both directions to ensure accuracy. Chromatograms were visually inspected and corrected as necessary in Geneious v. 11 (Biomatters). COI alignments were translated into amino acids to check for stop codons and frameshift mutations, then trimmed to 658 bp. 16S rRNA sequences were aligned with MAFFT v. 7 [38, 39] using the E-INS-i algorithm which performed better on our expanded dataset than the Q-INS-i algorithm used by Wiggering et al.  in aligning variable regions. All newly generated sequences have been deposited in GenBank. See Table 1 for sequence accession numbers and voucher information.
Species delimitation analyses
The final dataset for Planaxis sulcatus comprised 88 COI and 81 16S rRNA sequences. Primary Species Hypotheses (PSHs) (sensu [40, 41]) were formulated using ASAP  for the COI dataset (default parameters) with a distance matrix constructed in PAUP* v.4.0  using the best-fit nucleotide substitution model (GTR + I + G) calculated with jModelTest v. 2 . ASAP is the updated implementation of the ABGD hierarchical clustering algorithm  which has the added functionality of calculating a so-called “ASAP-score” that is used to rank the alternative partitioning schemes. It also has the benefit of obviating the need in ABGD for an a priori, user-defined range of P values which correspond to the minimum and maximum intraspecific divergence. Wiggering et al. defined Pmin = 0.0031 and Pmax = 0.041 based on COI distance values observed in a distantly related gastropod genus, Littoraria (Littorinoidea).
PSHs were then tested for reciprocal monophyly [45, 46] through phylogenetic analysis on the combined dataset (COI partitioned by codon position + 16S rRNA), using Bayesian (BA) and maximum likelihood (ML) methods. Because our dataset includes a substantial number of additional specimens compared with that of Wiggering et al., we re-estimated the best-fit nucleotide substitution models. Nucleotide substitution models were selected using jModelTest and were as follows: COI 1st codon position = SYM + I; COI 2nd codon position = F81; COI 3rd codon position = GTR + G; 16S rRNA = GTR + I + G. Bayesian analyses were run using MrBayes v. 3.2  on the CIPRES Science Gateway v. 3.3  (107 generation, 25% burn-in). MCMC convergence was assumed to be reached when the effective sample size was >200 and values of the potential scale reduction factor were 1 (analysed with Tracer v. 1.7 ). ML analyses were run online using W-IQ-TREE v. 1  (ultrafast bootstrap replicates = 1000). A node was considered supported if the corresponding Bayesian posterior probability (PP) was ≥0.95. The same threshold was used for interpreting the ultrafast bootstrap value (UFb) obtained from the ML analysis, as suggested by the authors of IQ-TREE .
PSHs were further explored with bGMYC , a Bayesian implementation of the general mixed Yule-coalescent model for species delimitation, performed on the COI dataset using the bGMYC package in R v. 3.2.1 . Ultrametric trees were generated using BEAST v.1.8  (2 runs of 208 generations, sampled each 1000, 25% burn-in, substitution model = HKY, clock model = lognormal relaxed clock, tree prior = Birth-Death Process). MCMC convergence was assessed with Tracer 1.7  and assumed to have occurred if effective sample sizes were greater than 200. The bGMYC analysis was run on 100 trees sampled equidistantly from those obtained from the BEAST analysis (generations = 50000, burn-in=40000, t1 = 1, t2 = 88, thinning = 100).
COI pairwise genetic distances were calculated using MEGA v.7  (nucleotide substitution model = Kimura 2-parameter, pairwise deletion of missing data). Median-joining  (MJ) haplotype networks were inferred and visualized using PopART (http://popart.otago.ac.nz).
The 10 best partitions found by ASAP divided the COI dataset into 2 to 56 hypothetical species; no partition included all specimens in a single species. The best ASAP partition divided the dataset into four PSHs, corresponding to Clades I, IIa, IIb+IIc, and III (Figs. 2, 3, and Additional file 1: Fig. S1; Table 1; clade names correspond to those used by Wiggering et al.). This was the best partition according to both ASAP P-rank and W score (5.22e−03 and 2.30e−03, respectively). The second-best partition according to the P-rank score divided the dataset into two PSHs (Clades I and II+III), while the second-best according to the W score identified five PSHs (Clades I, IIa, IIb, IIc, and III). All subsequent partitions identified between 25 and 56 PSHs. The bGMYC analysis (Fig. 2 and Additional file 1: Fig. S2) recovered five partitions that divided the dataset into two to six groups, with increasing posterior probability values. The partition with the highest PP value (PP=0.4) divided the COI dataset into six groups, considering specimens MNHN-IM-2019-16132, UF 4639999, and ZMB107725-6 as an entity distinct from Clade I. The other partitions divided the dataset into the same clades as found by ASAP
All phylogenetic analyses supported monophyly of Planaxis sulcatus sensu lato (Fig. 2, Additional file 1: Figs. S3–S6; PP=1 UFb=98). In analyses of the concatenated COI+16S rRNA dataset, monophyly of Clades II+III (PP=1 UFb=96), IIa (PP=1 Ufb=100), and IIc (PP=0.97 UFb=98) was supported by both BA and ML analyses. Monophyly of Clades I and IIb+IIc was supported by only one of the two methods (PP=1 and PP=0.95, respectively), while that of Clade IIb and III was not supported (PP=0.79 UFb=91, and PP=0.81 UFb=93, respectively). The BA tree resulting from analysis of the COI dataset recovered two reciprocally monophyletic clades (I and II+III, both PP=1), while all other single-gene analyses with BA or ML did not support any of the other subdivisions within P. sulcatus.
COI genetic distances within and between clades were calculated using K2P model-corrected distances (Table 2). For the two-PSH partition, the range of intra-clade genetic distances was 0–6.19%, and 6.32–12.09% between clades. For the four- and five-PSH partitions, the ranges of intra-clade genetic distances were 0–4.02% and 0–2.56% respectively, while the ranges of inter-clade distances were almost identical for both (3.69–12.09% and 3.68–12.09%, respectively). The range of genetic distances when all P. sulcatus specimens were considered to belong to a single species was 0–12.09%.
The MJ network analysis divided P. sulcatus haplotypes into two macro haplogroups corresponding to Clades I and II+III (Fig. 4), separated by at least 13 polymorphic sites. Haplotypes corresponding to Clades IIa, IIb, IIc, and III were each separated by between 8 to 9 polymorphic sites.
One, two, or five species?
Our goal is not to resolve with certainty the number of distinct evolutionary lineages that comprise the Planaxis sulcatus complex. This will require an integrative taxonomic approach including detailed morphological (shell, radula) comparisons of name-bearing types in historical collections with the sequenced specimens. However, it seems certain that the evidence supports the interpretation that P. sulcatus comprises at least two evolutionary lineages, and quite probably more than that. The existence of at least two distinct lineages comprising Clade I (Central Indo-Pacific, North-western Pacific, and Palmyra Atoll), and Clade II (Red Sea, Arabian Sea, Persian Gulf, and Western Indian Ocean) + Clade III (Western Indian Ocean, and Central Indo-Pacific) was supported by the second-best ASAP P-rank score, one of the bGMYC partitions, reciprocal monophyly in phylogenetic analyses of the COI and COI+16S rRNA datasets, and, the MJ network of COI sequences that identified two major distinct haplogroups.
In the Wiggering et al.  analysis, ABGD and bGMYC converged on the interpretation that three evolutionary lineages are present in their dataset, corresponding to Clades I, II and III. Here we have added sequences for 77 individuals, more than tripling the size of the original dataset, and have expanded the geographic coverage to include southern Madagascar, South Africa, Western Australia, Singapore, China, Japan, Papua New Guinea, Philippines, Vanuatu and Palmyra. Interestingly, this tripartite subdivision is neither considered as the best partition by any of the species delimitation analyses of the expanded molecular dataset, nor supported by the phylogenetic analyses (the monophyly of Clade III was not supported).
Wiggering et al.  reported a maximum K2P model-corrected divergence in COI of 3.08% and interpreted it as indicative of intraspecific population structuring. However, this value is the result of a misinterpretation of the ABGD methodology. ABGD uses a range of prior intraspecific divergences (P; also referred to as prior maximal distances) to partition a dataset into putative species based on a statistically inferred barcode gap. ABGD uses a value of P within a user-defined range, which Wiggering et al. specified as Pmin = 0.0031 and Pmax = 0.0411, and searches for a barcoding gap above this threshold. Thus, the value of 3.08% reported in Wiggering et al. (Table S4) is the prior maximal distance that yields a single-species partition, and not the maximum K2P intraspecific divergence that can be calculated for the dataset. In fact, we recalculated a maximum value of 11.6% for the Wiggering et al. dataset, which is more consistent with the value of 12.09% obtained here for the expanded dataset. Both of these upper values far exceed what is considered as indicative of maximum species-level divergence for other marine caenogastropods (~ 2.2%-4.6%) [27, 57,58,59]. Even when considering the two-PSH partition recognizing Clades I and II+III, the maximum intra-clade genetic divergence (6.19%) still exceeds this threshold. The bGMYC analysis suggested the presence of up to six hypothetical species; however, none of the other methods used for the species delimitation analysis supported this hypothesis, and it may simply reflect the known tendency of bGMYC to oversplit the molecular dataset [60,61,62]. Our analyses produced two additional hypotheses, even if less supported, with four or five PSHs, identified as Clades I, IIa, IIb+IIc, and III (Fig. 2). These hypotheses were the best- and second best-scoring partitions according to ASAP P-rank and W score and yielded maximum intra-clade genetic distances in COI of 4.02% and 2.56%, respectively, still relatively high, but more consistent with hypothesized species-level thresholds for marine caenogastropods. The MJ haplotype network produced haplogroups that are less differentiated but still distinct. Moreover, high levels of inter-clade sequence divergence are maintained even where members of different clades occur in sympatry. For example, MNHN-IM-2009-27091 (Clade IIa), MNHN-IM-2009-27101 and MNHN-IM-2009-27120 (Clade III), collected at the same site within 100 m of each other in south Madagascar, differ by 5.02% sequence divergence in COI; ZMB107725-4 (Clade I) and ZMB107933-1 (Clade III), also collected at the same site, differed by 8.93%. In Indonesia, ZMB106461-1 and ZMB106003-h1, also from Clades I and III respectively, while not syntopic, were collected from two sites only ~10 km apart on Muna Island and differed by 7.49%. Sympatry of specimens from genetically divergent clades, for which there is no concordance in the results of different delimitation methods, as in the present case, is generally considered a diagnostic criterion for species delimitation (e.g. ). Even if some of our PSHs were not uniformly corroborated by the phylogenetic analyses, with some of the clades statistically supported only by one of the two methods in analyses of the concatenated dataset and were not always reciprocally monophyletic in the single-gene analyses, still the levels of genetic differentiation in sympatry suggest the presence of more than two species. The lack of support for Clades II and IIb, in particular, is undoubtedly a consequence of the comparatively poor sampling and the quantity of missing data; only 13 individuals were sequenced in Clade II and its subclades, and only four individuals were sequenced for both markers.
The case for poecilogony
Wiggering et al. stated that, even if Planaxis sulcatus comprised a complex of cryptic species, their claim of poecilogony is upheld in Clade III, for which they documented individuals bearing both veligers and juveniles. However, even here the evidence seems inconclusive.
Clade III is distributed in the Indian Ocean and western Pacific and, in the analysis of Wiggering et al., included 12 sequenced specimens from Indonesia, Fiji, Thailand, Mozambique, Madagascar, and Mauritius. However, inferences of developmental mode were based on observations of 12 gravid individuals, none of which were sequenced, and only some of which can be considered unambiguous.
The fact that none were sequenced is problematic given that five of the observations were from Indonesia, an area of overlap between Clades I and III. With no apparent way to differentiate morphologically between adults of different clades, it is not clear how mode of development was assigned to the sequenced terminals. Indeed, all observations of developmental mode must be considered suspect in the region of overlap between the two clades (i.e., Indonesia, Thailand, Fiji). In the case of Indonesia, this issue is of little consequence given that all examined gravid females had at most late larvae in their brood pouch and all terminals were inferred to possess veligers. In the case of Fiji, the terminal was coded based on the dissection of one gravid female, however, the expanded molecular analysis indicates that clades I and III may overlap in this area as well. Thus, it is not clear to which clade this developmental mode should be assigned.
Some geographic areas were not represented by gravid females, chief among them being Thailand and Madagascar. In these cases, Wiggering et al. “assumed reproductive mode based on area of origin.” This led to inferences of developmental mode for individuals in their tree that were not empirically based, and in some cases conflicted with observations. For example, this coding strategy resulted in two individuals collected at the same site in Thailand (Clade I, ZMB107725-4; Clade III, ZMB107933-1) being assigned to different developmental modes, veligers and juveniles, respectively, despite the fact that no gravid females were examined from Thailand. For Mauritius, three terminals from Clade III were coded with two different developmental modes (ZMB117931-2, ZMB117937 as veligers, and ZMB117931-1, as juveniles), with observations of only a single gravid female bearing eggs/early larvae. For Mozambique, four gravid females were observed, three with larvae and one with juveniles, none of which were sequenced, but both sequenced terminals were inferred to possess juveniles. This inference was apparently strengthened by one observation of a gravid female from Tanzania ostensibly bearing juveniles <0.5 mm, also with no corresponding sequence. However, the specimen examined from Tanzania (ZMB 108265-15), listed in Wiggering et al. Table S3 as bearing large numbers (6980) of juveniles, is figured in their Fig. 1c as a late larva.
Observations of larvae in Clade III included one individual from Fiji, five from Indonesia, one from Mauritius, and three from Mozambique. The six observations from Fiji and Indonesia were obtained from unsequenced individuals collected in the area of overlap between Clades I and III, and potentially, in the case of Fiji, from a mislocalized specimen (see “Future research”, below). Thus, these observations should be considered unreliable. This leaves four observations from Mauritius and Mozambique. The adage, absence of evidence is not evidence of absence, pertains here. As acknowledged by Wiggering et al., an observation of larvae in the brood pouch is inconclusive with regard to mode of development given the possibility that the individual was sampled early in the brooding phase. For the two observations of juveniles, one of them is erroneous, leaving the entire case of poecilogony in Clade III to rest on a single observation from Mozambique which was extrapolated to five terminals on their tree, sometimes in direct contradiction to observations of gravid females from the same geographic area.
Thus, the coding strategy used by Wiggering et al.  to assign developmental mode to terminals in their tree extrapolated inconsistently from a few observations for mostly unsequenced individuals; this made their tree difficult to interpret at best, and misleading at worst. Their claim that “developmental modes do not entirely correlate with genetic clades” (Wiggering et al. : 5) does not carry much weight.
As outlined by Knott & McHugh , the first step in identifying reliable cases of poecilogony is to convincingly rule out potentially cryptic species. Resolving this issue with certainty necessitates advancing several lines of investigation. The first will require determining the number of distinct evolutionary lineages and their rank within an integrative taxonomic framework. It is possible that morphological characters may come to light that allow discrimination of the molecular clades. For example, two recently identified molecular lineages formerly recognized as Supplanaxis nucleus in the Caribbean were shown through careful scrutiny to be morphologically diagnosable in features of the shell and radula . Indeed, Bandel  observed differences in radular morphology of Planaxis sulcatus individuals from the Red Sea compared to those from South Africa , and others from New Caledonia , representing three different molecular clades in our tree. Bandel  noted that some of the differences were of so great a magnitude that they could be indicative of species-level differences.
A related line of enquiry will require refining geographical distributions and observations of developmental mode in the north-western Indian Ocean. Assuming that Clade III represents a distinct evolutionary lineage, as described above, the case for poecilogony made by Wiggering et al. in this clade rests on essentially just one observation of an individual from northern Mozambique bearing juveniles 0.6–1.0 mm in size (AMS 322969-8; Wiggering et al. Table S3) in a clade where all other observations were only of larvae. One potential explanation is that this represents an extension of Clade II outside the Red Sea and Persian Gulf. Barkati and Ahmed  reported ‘direct’ development in individuals from Karachi, Pakistan which may represent an extension of Clade II to the northern Arabian Sea. Another potential explanation is that all observations of larvae in the Indian Ocean were made on individuals early in the brooding phase. However, the situation in the north-western Indian Ocean is far more complex than the analysis of Wiggering et al. would suggest. Although the Red Sea population has been recognized as a distinct species by Dekker & Orlin  and Janssen et al. , it is not known how far the Red Sea form is distributed in the northern Indian Ocean and, further, not all Red Sea Planaxis sulcatus produce juveniles. As reported by Hulings  and Bandel , individuals from the Gulf of Aqaba release veligers. The presence of two developmental modes in the Red Sea may support the interpretation that two evolutionary lineages are represented here, as indicated by the five-PSH partition. Thus, the case for or against poecilogony in Planaxis sulcatus will depend on the resolution of clade structure and distribution of developmental mode in the north-western Indian Ocean.
It should be noted that we consider the record from Fiji reported in Wiggering et al. from Clade III (ZMB 106376 h2) to be problematic. Our analyses indicated that Clade I (UF 295967, UF 330842) occurs in Fiji, which is more consistent with the interpretation of a single widespread clade in the western Pacific. The presence of Clade III there would require a gap spanning the comparatively well sampled areas of South East Asia, northern Australia, New Caledonia, and Vanuatu. Although our geographical coverage is far from comprehensive, this seems unlikely even with present sampling effort. Thus, we consider their record from Fiji as dubious pending further study, possibly representing a mislocalized specimen or sequencing contamination.
Lastly, Wiggering et al. claimed there are no morphological differences among larvae in the brood pouch, regardless of inferred development type. It may be difficult to assess such differences through light microscopic analysis, but it should be possible to recognize differences in mode of larval nourishment from the size of the embryonic shells. This likely requires higher magnification and higher resolution (e.g., scanning electron microscopy) to assess. Moreover, while differences in the embryonic shells may be subtle, mature veliger larvae that have been brooded (Wiggering et al. : Fig. 1d, E) are readily differentiable from free-swimming forms in size, ornament and number of whorls; differences in size and number of whorls are also evident among free-swimming veligers from different clades (see Houbrick : fig. 1e; Bandel : pl. 2, figs: 5, 7, 9). Unfortunately, adult planaxids do not retain undamaged protoconchs into adulthood given the unforgiving habitats where they live, and postlarval juveniles are rare in collections. Thus, documenting protoconch morphology from specimens in historical collections is a challenge. Future collecting strategies should include subadults and juveniles as possible.
We conclude that the evidence for the existence of a single, widespread, species of Planaxis sulcatus, and for the presence of poecilogony in Clade III, is equivocal. The hypothesis that Planaxis sulcatus comprises at least two distinct evolutionary lineages was instead robust to the addition of a significant number of new sequences and to expanding the geographic coverage. Additional sampling of specimens is needed, especially in Clade II, as well as a broader comparative analysis of the protoconchs to assess morphological differences within and among the clades, particularly of the mature veliger stages, as part of an integrative taxonomic approach. Taxonomy of this complex would certainly benefit from the use of additional molecular markers, including nuclear ones (e.g., ribosomal intergenic spacers, microsatellites, RAD markers) to reach a more robust species delimitation.
Availability of data and materials
The genetic sequences generated during the current study are available in the GenBank repository (MZ470478–MZ470583). Permanent link to sequences: https://www.ncbi.nlm.nih.gov/nuccore/?term=MZ470528:MZ470583%5Baccn%5D
Barcode of Life Data System
Muséum national d’Histoire naturelle
National Museum of Natural History
Florida Museum of Natural History
Museum für Naturkunde
Haszprunar G, Salvini-Plawen L, Rieger RM. Larval planktotrophy—a primitive trait in the Bilateria? Acta Zool. 1995;76:141–54.
Strathmann RR. Feeding and nonfeeding larval development and life-history evolution in marine invertebrates. Annu Rev Ecol Syst. 1985;16:339–61.
Jablonski D, Lutz RA. Larval ecology of marine benthic invertebrates: paleobiological implications. Biol Rev. 1983;58:21–89. https://doi.org/10.1111/j.1469-185X.1983.tb00380.x.
Shuto T. Larval ecology of prosobranch gastropods and its bearing on biogeography and paleontology. Lethaia. 1974;7:239–56.
Oliverio M. Life-histories, speciation, and biodiversity in Mediterranean prosobranch gastropods. Vie Milieu. 1996;46:163–9.
Giard AM. La poecilogonie. Comp Rend Six Congr Int Zool Berne. 1904;1905:617–46.
Chia F-S, Gibson G, Qian P-Y. Poecilogony as a reproductive strategy of marine invertebrates. Oceanol Acta. 1996;19:203–8.
Krug PJ. Poecilogony and larval ecology in the gastropod genus Alderia. Am Malacol Bull. 2007;23:99–111.
Krug PJ. Not my “type”: larval dispersal dimorphisms and bet-hedging in opisthobranch life histories. Biol Bull. 2009;216:355–72.
Ellingson RA, Krug PJ. Evolution of poecilogony from planktotrophy: cryptic speciation, phylogeography, and larval development in the gastropod genus Alderia. Evolution (N Y). 2006;60:2293–310.
Vendetti JE, Trowbridge CD, Krug PJ. Poecilogony and population genetic structure in Elysia pusilla (Heterobranchia: Sacoglossa), and reproductive data for five sacoglossans that express dimorphisms in larval development. Integr Comp Biol. 2012;52:138–50.
Blake JA, Arnofsky PL. Reproduction and larval development of the spioniform Polychaeta with application to systematics and phylogeny. Hydrobiologia. 1999;402:57–106.
Hoagland KE, Robertson R. An assessment of poecilogony in marine invertebrates: phenomenon or fantasy? Biol Bull. 1988;174:109–25.
Bouchet P. A review of poecilogony in gastropods. J Molluscan Stud. 1989;55:67–78.
Mcdonald KA, Collin R, Lesoway MP. Poecilogony in the caenogastropod Calyptraea lichen (Mollusca: Gastropoda). Invertebr Biol. 2014;133:213–20.
Russini V, Giannuzzi-Savelli R, Pusateri F, Prkić J, Fassio G, Modica MV, et al. Candidate cases of poecilogony in Neogastropoda: implications for the systematics of the genus Raphitoma Bellardi, 1847. Invertebr Syst. 2020;34:293–318.
Risbec J. Biologie et ponte de Mollusques gastéropodes néo-calédoniens. Bull Société Zool Fr. 1935;60:387–417.
Thorson G. Studies on the egg masses and larval development of Gastropoda from the Iranian Gulf. Danish Sci Investig Iran. 1940;2:159–238.
Wiggering B, Neiber MT, Gebauer K, Glaubrecht M. One species, two developmental modes: a case of geographic poecilogony in marine gastropods. BMC Evol Biol. 2020;20:76.
Knott KE, McHugh D. Introduction to Symposium: Poecilogony–a window on larval evolutionary transitions in marine invertebrates. Integr Comp Biol. 2012;52:120–7.
Spalding MD, Fox HE, Allen GR, Davidson N, Ferdaña ZA, Finlayson M, et al. Marine ecoregions of the world: a bioregionalization of coastal and shelf areas. Bioscience. 2007;57:573–83.
Lydeard C, Holznagel WE, Glaubrecht M, Ponder WF. Molecular phylogeny of a circum-global, diverse gastropod superfamily (Cerithioidea: Mollusca: Caenogastropoda): pushing the deepest phylogenetic limits of mitochondrial LSU rDNA sequences. Mol Phylogenet Evol. 2002;22:399–406.
Ran K, Li Q, Qi L, Li W, Kong L. DNA barcoding for identification of marine gastropod species from Hainan island China. Fish Res. 2020;225:105504.
Zou S, Li Q, Kong L. Additional gene data and increased sampling give new insights into the phylogenetic relationships of Neogastropoda, within the caenogastropod phylogenetic framework. Mol Phylogenet Evol. 2011;61:425–35.
Setiamarga DHE, Nakaji N, Iwamoto S, Teruya S, Sasaki T. DNA barcoding study of shelled gastropods In the intertidal rocky coasts of central Wakayama Prefecture, Japan, using two gene markers. Int J GEOMATE. 2019;17:9–16.
Ozawa T, Kohler F, Reid DG, Glaubrecht M. Tethyan relicts on continental coastlines of the northwestern Pacific Ocean and Australasia: molecular phylogeny and fossil record of batillariid gastropods (Caenogastropoda, Cerithioidea). Zool Scr. 2009;38:503–25.
Sun Y, Li Q, Kong L, Zheng X. DNA barcoding of Caenogastropoda along coast of China based on the COI gene. Mol Ecol Resour. 2012;12:209–18.
Izadian M. The molecular identification of four species of Gastropoda on rocky shores of the Persian Gulf. Islam Azad Univ. 2020;11:59–74.
Ardura A, Planes S, Garcia-Vazquez E. Aliens in paradise. Boat density and exotic coastal mollusks in Moorea Island (French Polynesia). Mar Environ Res. 2015;112:56–63. https://doi.org/10.1016/j.marenvres.2015.08.007.
Galindo LA, Puillandre N, Strong EE, Bouchet P. Using microwaves to prepare gastropods for DNA barcoding. Mol Ecol Resour. 2014;14:700–5.
Geller J, Meyer C, Parker M, Hawk H. Redesign of PCR primers for mitochondrial cytochrome c oxidase subunit I for marine invertebrates and application in all-taxa biotic surveys. Mol Ecol Resour. 2013;13:851–61.
Strong EE, Whelan NV. Assessing the diversity of western North American Juga (Semisulcospiridae, Gastropoda). Mol Phylogenet Evol. 2019;136:87–103.
Palumbi S, Martin A, Romano S, McMillan W, L S, Grabowski G. The simple fool’s guide to PCR, version 2.0. privately published by S. Palumbi, Department of Zoology, University of Hawaii, Honolulu, HI; 1991.
Aljanabi SM, Martinez I. Universal and rapid salt-extraction of high quality genomic DNA for PCR-based techniques. Nucleic Acids Res. 1997;25(22):4692–3.
Palumbi SR. Nucleic acids II: the polymerase chain reaction. Mol Syst. 1996;2:205–47.
Hayashi S. The molecular phylogeny of the Buccinidae (Caenogastropoda: Neogastropoda) as inferred from the complete mitochondrial 16S rRNA gene sequences of selected representatives. Molluscan Res. 2003;25:85–98.
Leray M, Yang JY, Meyer CP, Mills SC, Agudelo N, Ranwez V, et al. A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: application for characterizing coral reef fish gut contents. Front Zool. 2013;10:34. https://doi.org/10.1186/1742-9994-10-34.
Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30:772–80.
Katoh K, Rozewicki J, Yamada KD. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 2019;20:1160–6. https://doi.org/10.1093/bib/bbx108.
Puillandre N, Modica MV, Zhang Y, Sirovich L, Boisselier M-C, Cruaud C, et al. Large-scale species delimitation method for hyperdiverse groups. Mol Ecol. 2012;21:2671–91. https://doi.org/10.1111/j.1365-294X.2012.05559.x.
Puillandre N, Lambert A, Brouillet S, Achaz G. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Mol Ecol. 2012;21:1864–77.
Puillandre N, Brouillet S, Achaz G. ASAP: assemble species by automatic partitioning. Mol Ecol Resour. 2020;21:609–20. https://doi.org/10.1111/1755-0998.13281.
Swofford DL. PAUP*. Phylogenetic Analysis Using Parsimony *and other methods. Version 4. Sunderland, Massachusetts: Sinauer Associates; 2002.
Darriba D, Posada D. jModelTest 2.0 Manual. Nat Methods. 2012;9:772. http://code.google.com/p/jmodeltest2.
Knowlton N. Molecular genetic analyses of species boundaries in the sea. Hydrobiologia. 2000;420:73–90.
Padula V, Bahia J, Stöger I, Camacho-García Y, Malaquias MAE, Cervera JL, et al. A test of color-based taxonomy in nudibranchs: Molecular phylogeny and species delimitation of the Felimida clenchi (Mollusca: Chromodorididae) species complex. Mol Phylogenet Evol. 2016;103:215–29. https://doi.org/10.1016/j.ympev.2016.07.019.
Ronquist F, Teslenko M, Van Der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61:539–42.
Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. In: Proceedings of the Gateway Computing Environments Workshop (GCE). New Orleans, LA; 2010. p. 1–8. doi: https://doi.org/10.1109/GCE.2010.5676129.
Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst Biol. 2018;67:901–4.
Trifinopoulos J, Nguyen L-T, von Haeseler A, Minh BQ. W-IQ-TREE: a fast online phylogenetic tool for maximum likelihood analysis. Nucleic Acids Res. 2016;44:W232–5. https://doi.org/10.1093/nar/gkw256.
Minh BQ, Nguyen MAT, von Haeseler A. Ultrafast approximation for phylogenetic bootstrap. Mol Biol Evol. 2013;30:1188–95. https://doi.org/10.1093/molbev/mst024.
Reid NM, Carstens BC. Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yule-coalescent model. BMC Evol Biol. 2012;12:196.
Team RC. R: A language and environment for statistical computing. 2014. http://www.r-project.org/.
Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214. https://doi.org/10.1186/1471-2148-7-214.
Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:msw054. https://doi.org/10.1093/molbev/msw054.
Bandelt HJ, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.
Collin R. Development, phylogeny, and taxonomy of Bostrycapulus (Caenogastropoda: Calyptraeidae), an ancient cryptic radiation. Zool J Linn Soc. 2005;144:75–101.
Zou S, Li Q, Kong L, Yu H, Zheng X. Comparing the usefulness of distance, monophyly and character-based DNA barcoding methods in species identification: a case study of Neogastropoda. PLoS ONE. 2011;6:e26619.
Russini V, Fassio G, Modica M, DeMaintenon M, Oliverio M. An assessment of the genus Columbella Lamarck, 1799 (Gastropoda: Columbellidae) from eastern Atlantic. Zoosystema. 2017;39:197–212.
Paz A, Crawford AJ. Molecular-based rapid inventories of sympatric diversity: a comparison of DNA barcode clustering methods applied to geography- based vs clade-based sampling of amphibians. J Biosci. 2012;37:887–96.
Pentinsaari M, Vos R, Mutanen M. Algorithmic single-locus species delimitation: effects of sampling effort, variation and nonmonophyly in four methods and 1870 species of beetles. Mol Ecol Resour. 2017;17:393–404.
Lin X-L, Stur E, Ekrem T. Exploring species boundaries with multiple genetic loci using empirical data from non- biting midges. Zool Scr. 2018;47:325–41.
Kekkonen M, Hebert PDN. DNA barcode-based delineation of putative species: efficient start for taxonomic workflows. Mol Ecol Resour. 2014;14:706–15.
Strong EE, Bouchet P. Hidden in plain sight: two co-occurring cryptic species of Supplanaxis in the Caribbean (Cerithioidea, Planaxidae). Zookeys. 2020;991:85.
Bandel K. The radulae of Caribbean and other Mesogastropoda and Neogastropoda. Zool Verh. 1984;214:1–199.
Barnard KH. Contributions to the knowledge of South African Mollusca III. Gastropoda: Prosobranchia: Taenioglossa. Ann South African Mus. 1963;47:1–199.
Barkati S, Ahmed M. Studies on the reproductive biology of some prosobranchs from the coast of Karachi (Pakistan) bordering the northern Arabian Sea, 1: Observations on Planaxis sulcatus (Born, 1780). The Veliger. 1982;24:355–8.
Dekker H, Orlin Z. Check-list of Red Sea Mollusca. Spirula. 2000;47(Suppl):3–46.
Janssen R, Zuschin M, Baal C. Gastropods and their habitats from the northern Red Sea (Egypt: Safaga). Part 2: Caenogastropoda: Sorbeoconcha and Littorinimorpha. Ann Naturhist Mus Wien, Ser A. 2011;113:373–509.
Hulings NC. Aspects of the reproduction of rocky intertidal mollusks from the Jordan Gulf of Aqaba (Red Sea). The Veliger. 1986;28:318–27.
Bandel K. Families of the Cerithioidea and related superfamilies (Palaeo-Caenogastropoda; Mollusca) from the Triassic to the Recent characterized by protoconch morphology - including the description of new taxa. Freib Forschungshefte C. 2006;511:59–138.
Houbrick RS. Anatomy, reproductive biology, and phylogeny of the Planaxidae (Cerithiacea: Prosobranchia). Washington, D.C.: Smithsonian Institution Press; 1987; Smithsonian Contributions to Zoology, 445: 1–57.
Many of the specimens sequenced herein were collected under the auspices of the Our Planet Reviewed program at the MNHN. We extend our thanks to John Slapcinsky (Florida Museum of Natural History) for providing six individuals preserved for molecular analysis, Benoit Dayrat (Pennsylvania State University) for kindly providing previously unpublished COI sequences from three individuals, and Dai Herbert and Herman van der Bank for unpublished COI sequences from four individuals from the South African barcoding project. We are grateful to the staff of Laboratories of Analytical Biology and Department of Invertebrate Zoology for generating the sequences and curating the specimens produced in Washington DC. Thanks also to Barbara Buge (MNHN) for curating the barcoding vouchers in Paris and to Nicolas Puillandre (MNHN) for discussions on the ABGD method. Two anonymous reviewers provided very useful comments that helped us in preparing a clearer revised version.
Ethics approval and consent to participate
Consent for publication
All authors gave consent for publication of the present work.
All authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Fassio, G., Bouchet, P., Oliverio, M. et al. Re-evaluating the case for poecilogony in the gastropod Planaxis sulcatus (Cerithioidea, Planaxidae). BMC Ecol Evo 22, 13 (2022). https://doi.org/10.1186/s12862-022-01961-7
- Reproductive biology
- Cryptic species
- Larval development