A reliable taxonomy, together with more accurate knowledge of the geographical distribution of species, is a fundamental element for the study of biodiversity. Multiple studies on the gastropod family Neritidae record three species of the genus Neritina in the Brazilian Province: Neritina zebra (Bruguière, 1792), Neritina virginea (Linnaeus, 1758), and Neritina meleagris Lamarck, 1822. While N. zebra has a well-established taxonomic status and geographical distribution, the same cannot be said regarding its congeners. A widely cited reference for the group in Brazil considers N. meleagris a junior synonym of N. virginea. Using a molecular approach (phylogenetic, species delimitation, and statistical parsimony network analyses), based on two mitochondrial markers (COI and 16S), this study investigated if N. virginea and N. meleagris are distinct species. The molecular results confirmed the existence of two strongly supported distinct taxonomic entities in the Brazilian Province, which is consistent with the morphological descriptions previously proposed for N. virginea and N. meleagris. These species occur in sympatry in the intertidal sandstone formations of Northeastern Brazil. Despite the great variation in the colour patterns of the shells, the present study reinforced previous observations that allowed the differentiation of these two species based on these patterns. It also emphasized the importance of the separation of these two clades in future studies, especially those conducted in the Brazilian Province, since these species may cohabit.

Brazilian Province, Caribbean Province, geographic distribution, neritids, species delimitation

Molluscs from the gastropod family Neritidae are the most diverse members of Neritimorpha (Kano et al. 2002), with some groups within this family having variable shell colouration patterns (e.g., Russell 1941; Tan and Clements 2008; Eichhorst 2016). Due to the great variety of colour patterns, the delimitation of different species could be hampered, especially if they are closely related and live in sympatry (e.g., Huang 1997; Blanco et al. 2014). This may explain the disparate estimates in the literature of the number of species of Neritina reported for the Brazilian Province.

Several studies report three species of the genus Neritina on the Brazilian coast: Neritina zebra (Bruguière, 1792), Neritina virginea (Linnaeus, 1758), and Neritina meleagris Lamarck, 1822 (e.g., Baker 1923; Russell 1941; Rios 1975; Matthews-Cascon et al. 1990; Díaz and Puyana 1994; Quintero-Galvis and Castro 2013; Eichhorst 2016). While N. zebra has a well-established taxonomic status and geographical distribution (Matthews-Cascon et al. 1990; Rios 2009; Barroso et al. 2012; Eichhorst 2016), there is uncertainty regarding its congeners. The shell catalogues of Rios (1985, 1994, 2009), a widely cited reference in studies conducted in Brazil, state that only two species occur in the Brazilian Province: N. virginea and N. zebra. In these compendia, N. meleagris is considered a junior synonym of N. virginea without any justification. Quintero-Galvis and Castro (2013), using a molecular phylogenetic approach to analyse specimens from the Colombian coast (Caribbean Province), concluded that N. meleagris and N. virginea are phylogenetically close, but different species. Since these species have a wide geographic distribution, encompassing the Caribbean and Brazilian Provinces (Barroso et al. 2016; Eichhorst 2016), that are separated by a recognized biogeographic barrier (the Amazon-Orinoco outflow) (Floeter et al. 2008; Barroso et al. 2016), the inclusion of specimens from both biogeographical provinces in phylogenetic analyses is desirable.

Since a reliable taxonomy, together with a more accurate knowledge about the geographical distribution of species, is fundamental to the study of biodiversity (Wheeler et al. 2004), the present study aims to investigate if N. virginea and N. meleagris are two distinct species, using molecular data (phylogenetic, species delimitation, and statistical parsimony network analyses).

We collected specimens from Barra Grande beach (Piauí State) (2°54.125'S, 41°24.573'W) and Camocim beach (Ceará State) (02°51.778'S, 41°51.57'W), both located in Northeastern Brazil, and preserved them in 95% ethanol. We identified the species using the literature (Russell 1941; Matthews-Cascon et al. 1990; Eichhorst 2016), primarily based on the shape and colour patterns of the shells. Specimens were collected under SISBIO permit no. 57473-3 and deposited in the malacological collection “Prof. Henry Ramos Matthews” – series B (CMPHRM-B) of Universidade Federal do Ceará (UFC). A total of 17 specimens, eight newly sequenced and nine already published by Quintero-Galvis and Castro (2013) and available on GenBank, were used for phylogenetic reconstruction (Table 1). All sequences used were attributed to nominal species considered valid in the literature (see Aktipis and Giribet 2010; Cook et al. 2010; Page et al. 2013; Quintero-Galvis and Castro 2013).

We extracted whole genomic DNA from the foot muscle of specimens, using the Qiagen DNeasy Blood & Tissue Kit (Qiagen, Valencia, CA, USA). The quality and integrity of the DNA obtained were evaluated in a micro-volume spectrophotometer. Amplification of double-stranded fragments from the cytochrome c oxidase I (COI) and 16S mitochondrial genes was achieved by polymerase chain reaction (PCR) using newly developed neritid-specific custom primers for the 16S gene [(16SNer_F 5’ACTACTCCGCCTGTTTATCAAA3’) and (16SNer_R 5’GGGCTTAAACCTAATGCACTT3’)] and modified versions of Folmer et al. (1994) primers for the COI gene [(LCO1490_mod 5’ATTCTACGAATCAYAAAGAYATTGG3’) and (HCO2198_mod 5’TAWACTTCAGGATGACCRAAAAATCA3’)]. The PCR was carried out using GoTaq Green Master Mix (Promega Corporation), 1.25 μL of each primer (10 μM stock), and 100 ng of DNA template in a 25 μL reaction volume. The PCR cycles for COI and 16S amplification consisted of an initial denaturation step at 95 °C for 2 min, followed by 35 cycles of denaturation at 95 °C for 30 s, annealing at 48–49 °C for 45 s and extension at 72 °C for 1 min, and a final extension at 72 °C for 5 min. The PCR products were then examined using gel electrophoresis on 1.3% Tris-Borate-EDTA-agarose gel stained with SYBR Safe DNA Gel Stain (Invitrogen). The PCR products showing strong bands in gel electrophoresis were purified with IllustraExoProStar - 1 Step (GE Healthcare Life Sciences), following its standard protocol, and sent for Sanger sequencing (Macrogen Inc., South Korea).

The forward and reverse sequences for each gene fragment were edited using Geneious v. 7.1 (Biomatters). The concatenated alignments of COI and 16S were conducted using the MAFFT program with the G-INS-I algorithm (Katoh and Standley 2013) using the default parameters, with additional inspection by eye for accuracy (see Suppl. material 1). As our 16S sequences (650 bp for N. virginea and 651 bp for N. meleagris) were longer than those available on GenBank, we used a minor homologous region in the alignment of this gene. However, we deposited the full 650–651 bp 16S sequences in GenBank. The combined dataset contained 1124 bp (639 bp for COI and 485 bp for 16S). Evolutionary relationships were estimated for the concatenated genetic markers using Bayesian inference (BI) and maximum likelihood (ML) analyses. The best-fit evolution models were determined using PartitionFinder (Lanfear et al. 2012), considering the positions of the codon for the COI gene, which codes for protein, and a single partition for the 16S gene. The corrected Bayesian Information Criterion (BIC) was used to select among the options. PartitionFinder selected respectively GTR + I, F81, and HKY + G as the best model for the three positions of the codon in COI, and GTR + G for 16S. Bayesian inference, using the previously mentioned partitions and models, was performed using the MrBayes program (Ronquist et al. 2012) and the dataset was run for 3 × 10⁷ generations, with Markov chains sampled every 1000 generations, and the standard 25% burn-in calculated. Convergence was checked using Tracer 1.6 (http://beast.bio.ed.ac.uk/Tracer). Tree branches were considered strongly supported if posterior probabilities were ≥ 0.90. Randomized accelerated maximum likelihood (RAxML) (Stamatakis 2006) was used to generate a ML tree with partitions under the evolution model GTR + G and with 1 × 10⁴ replications. Branches with bootstrap values greater than 70 were considered strongly supported. Phylogenetic trees were drawn and edited in FigTree 1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

For the species delimitation analyses, we initially constructed a distance matrix based on the Kimura 2-parameter (K2P) model, using the COI sequences, in the MEGA 6.0.6 software (Tamura et al. 2013). This matrix was analysed with the default settings of the Automatic Barcode Gap Discovery method (ABGD) (Puillandre et al. 2012) (available at http://www.abi.snv.jussieu.fr/public/abgd/abgdweb.html). We also used the Species Delimitation plugin v1.04 for Geneious v. 7.1 (Masters et al. 2011) with two data sets: (1) the results of our Bayesian concatenated phylogenetic analysis (COI + 16S), and (2) the results of a neighbor-joining tree based on the K2P model with 10,000 bootstrap replicates using COI sequences generated in MEGA 6.0.6. In this analysis, we calculated (1) the mean distance between the members within the clade (Intra Dist), (2) the mean distance of those individuals to the nearest clade (Inter Dist-closest), (3) the ratio between Intra Dist and Inter Dist-closest, and (4) the P ID, which represents the mean probability (95% confidence interval) of correctly identifying an unknown member of the putative species to fit inside (Strict P ID), or at least to be the sister group of (Liberal P ID), the species clade in a tree (Masters et al. 2011).

A statistical parsimony network analysis was conducted with COI sequences (347 bp), using the TCS algorithm (Clement et al. 2002) implemented in PopART v. 1.7.2 (Leigh and Bryant 2015). The sequence alignment step followed the same procedures already described in our phylogenetic analysis protocol (see Suppl. material 1). In addition to the N. virginea and N. meleagris sequences generated in this study, we included 55 sequences of N. virginea from island (Puerto Rico: 44 sequences) and continental (Panama: 10; Colombia: 1) locations in the Caribbean Province (Aktipis and Giribet 2010; Cook et al. 2010; Page et al. 2013; Quintero-Galvis and Castro 2013) (Table 1). Along with phylogenetic and species delimitation analyses, we also included the only COI sequence available on GenBank attributed to N. meleagris from Colombia (Quintero-Galvis and Castro 2013) (Table 1).

The shells and opercula of the specimens of N. virginea and N. meleagris submitted to molecular analyses were observed and photographed under a stereomicroscope. A scanning electron microscope (SEM) was used to view their radulae (two females of N. virginea and two males and one female of N. meleagris) in the Analytical Facility of UFC (Central Analítica, UFC). This information was collected in order to compare our results with the information available in the literature.

The results of our molecular analyses (phylogeny, species delimitation, and statistical parsimony network) confirmed the existence of two strongly supported clades living in sympatry in the intertidal beachrocks of Northeastern Brazil (Brazilian Province). The Bayesian and maximum likelihood trees showed the same topology, with the formation of four clades within Neritina: Group I (Neritina piratica + Neritina usnea), Group II (Neritina meleagris, collected in NE Brazil), Group III (Neritina virginea, collected in NE Brazil and Colombia, + “Neritina meleagris”, from Colombia), and Group IV (Neritina punctulata) (Fig. 1). The Groups II and III were also observed in the ABGD analysis. Intraspecific distances in these groups were at least one order of magnitude smaller than the interspecific distances. Our minimum interspecific genetic distance values (COI region only), involving groups II and III, were 8.4 and 9.6%, respectively (Table 2). These values are higher than the minimum value assumed by Abdou et al. (2017) to characterize distinct Indo-Pacific Neritina species. In addition, the probabilities of a new sequence fitting inside P ID (Strict) or at least the sister group P ID (Liberal) of these clades were equal to, or in most cases, greater than 84% (Table 2). These results are compatible with the values found to delimit species in different groups of gastropods (e.g. Churchill et al. 2014; Cooke et al. 2014; Espinoza et al. 2014).

Figure 1. 

Molecular phylogenetic hypothesis (Bayesian tree) of some species of Neritidae of the Western Atlantic. The Bayesian tree was based on partial mitochondrial COI and 16S sequences. The Neritina meleagris and Neritina virginea clades (ingroup) are highlighted. The other taxa were used as outgroup. Numbers on and below the main branches represent the posterior Bayesian probabilities (BP) (>0.90) and bootstrap values for maximum likelihood (ML) (>70%), respectively. Specimens with the number “1” are from Camocim beach (Ceará State, NE Brazil) and those with numbers “2”, “3”, and “4” are from Barra Grande beach (Piauí State, NE Brazil). The numbered specimens of N. virginea (1, 2, 3, and 4) and N. meleagris (1, 2, 3, and 4) are the same specimens shown in Figure 3.

Although the distinction between clades showed high support values, the phylogenetic relationship between them could not be recovered. As we did not have access to the specimens, it was not possible to check the shell colour patterns of the Neritina meleagris from Colombia (obtained from GenBank) that was included in the same clade as Neritina virginea. Thus, we suspect that an error may have occurred at the time of submission of the sequences to GenBank, since, in the study of Quintero-Galvis and Castro (2013), N. virginea and N. meleagris appeared in very distinct branches of the phylogenetic tree.

Despite the geographical distance, all N. virginea sequences from Brazil and the Caribbean were very similar, with all haplotypes grouped within a few mutational steps (Fig. 2). This result reinforces the validity of N. virginea and confirms its presence in the Brazilian Province. As also observed in the phylogenetic analysis, the only sequence assigned as N. meleagris from the Caribbean Province is positioned within one of the most frequent N. virginea haplotypes for this region (Fig. 2). With respect to our N. meleagris sequences, although this species is found in sympatry with its congener in the Brazilian Province, it is separated from the N. virginea haplogroups by at least 36 mutational steps.

Figure 2. 

Statistical parsimony network analysis (TCS algorithm) based on 64 partial mitochondrial COI sequences (347 bp). This analysis included specimens of Neritina meleagris and Neritina virginea from the Caribbean and Brazilian Provinces. Size of the circle is proportional to frequency of the haplotype and colours inside the circles designate geographical locations to which the samples belong. Black circles correspond to hypothetical haplotypes. The number of mutational steps is indicated by dashes on branches. We highlighted the 36 mutational steps that separate the two species haplotypes.

Figure 3 shows the colour patterns, opercula, and radulae of the N. meleagris and N. virginea specimens collected from the Barra Grande and Camocim beaches. Our molecular results are consistent with the morphological descriptions previously proposed for each species (Russell 1941; Matthews-Cascon et al. 1990; Eichhorst 2016). Russell (1941) described, for both species, a colour pattern consisting of dark zigzag lines and lighter spots. However, this author emphasized that while N. virginea has a leading edge outlined in heavy black, N. meleagris instead has a leading edge outlined with white, white and black, or white and red, resembling imbricating scales. The imbricating scales pattern was emphasized in the original description of N. meleagris (Lamarck 1822), whereas Matthews-Cascon et al. (1990) and Eichhorst (2016) highlighted the differences in the leading edges of the colour pattern for each species. Although did not examine the type specimens, individuals of N. virginea from the Linnean Collection at the Natural History Museum, London (see http://linnean-online.org/17152/), and N. meleagris, from the type locality (Dominican Republic) (see http://data.biodiversitydata.nl/naturalis/specimen/ZMA.MOLL.313038), had the same leading edge patterns as described earlier. All analysed specimens of N. meleagris had the leading edge outlined with white or white and black, while N. virginea specimens had the leading edge outlined in black (Fig. 3A–H). Despite the great variation in their shell colour patterns, a more detailed observation of the leading edges of the N. virginea and N. meleagris shells allows the separation of the two species, even in the field. Warmke and Abbott (1962) also emphasized the use of leading edges to separate the two species. Williams (2017) argued that the colours and patterns of gastropod shells could be genetically determined, influenced directly by environmental factors, or a combination of both. Specifically, the patterns of leading edges (leading edge outlined with white or white and black in N. meleagris and outlined in black in N. virginea) appear to be under genetic control rather than be influenced directly by environmental factors, since the patterns for each species are consistent regardless of the location studied (e.g. Russell 1941; Eichhorst 2016; present study). This observation is reinforced by the clades obtained in our phylogenetic analysis, corroborating the diagnostic colour patterns previously described (Fig. 1).

Figure 3. 

Colour patterns of shells, opercula, and radulae of the Neritina virginea and Neritina meleagris analysed. The red arrows highlight the differences between the leading edges of colour patterns of both species: N. virginea has the leading edges outlined in heavy black, while N. meleagris has the leading edge outlined in white or black and white. A Neritina virginea_1 B Neritina virginea_2 C Neritina virginea_3 D Neritina virginea_4 E Neritina meleagris_1 F Neritina meleagris_2 G Neritina meleagris_3 H Neritina meleagris_4 I ventral view of shell of Neritina virginea J operculum (outer and inner views) of Neritina virginea K ventral view of shell of Neritina meleagris L operculum (outer and inner views) of Neritina meleagris M radula of Neritina virginea (SEM), with rachidian tooth enlarged in the upper left quadrant N radula of Neritina meleagris (SEM), with rachidian tooth enlarged in the upper left quadrant. Abbreviations: l1 first lateral tooth, l4 fourth lateral tooth, m marginal teeth, r rachidian tooth. The specimens with the number “1” are from Camocim beach (Ceará State, NE Brazil) and those with numbers “2”, “3”, and “4” are from Barra Grande beach (Piauí State, NE Brazil). The numbered specimens of N. virginea (1, 2, 3, and 4) and N. meleagris (1, 2, 3, and 4) are the same specimens used in the phylogenetic analysis of Figure 1. Scale bars: 1.0 mm (A–L); 100 μm (M, N).

Besides the shell colour patterns, N. virginea and N. meleagris differ from each other in subtle ways. The inner lips of the shells of the two species are denticulated. However, in N. virginea there are several small denticles interspersed by two larger teeth, while in N. meleagris the teeth are larger, more prominent in the central region, and less numerous when compared to N. virginea (Fig. 3I, K). Russell (1941), Matthews-Cascon et al. (1990), and Eichhorst (2016) also highlighted these differences regarding the number of teeth on the inner lip. Both species have a calcareous and smooth operculum, with a bifurcated apophysis. Comparing the opercula, N. virginea has a darker (bluish-black) and more elongated operculum, with the apophysis elements thinner and more separated from each other. On the other hand, the operculum of N. meleagris presents a lighter coloration (yellowish-black) and a semi-circular shape, with the apophysis elements much stouter and closer to each other (Fig. 3J, L). In the present study, Neritina virginea and N. meleagris have a very similar morphology of the radula: a rhipidoglossate radula, with one rachidian tooth, five pairs of lateral teeth, and many denticulated marginal teeth arranged in transverse rows (Fig. 3M, N; see also Suppl. material 2). The most striking difference between these radulae is the rectangular rachidian tooth, which has three cusps in N. meleagris (both male and female) but is cuspless in N. virginea. The first lateral tooth of N. virginea is more slender than that of N. meleagris. Previous studies have shown that the radula teeth pattern of neritids is very stable, the most variable character being the number of cusps on the fifth lateral tooth, which is likely correlated with age (Baker 1923; Huang 1997; Haynes 2001). This characteristic makes it difficult to define intra- and interspecific differences. Further studies are needed to better define the differences between the radulae of the two species.

Our molecular data show that N. virginea and N. meleagris are two distinct species, thus confirming the N. meleagris record for the Brazilian coast. In summary, our results, along with the already well-established record of Neritina zebra (Matthews-Cascon et al. 1990; Rios 2009; Barroso et al. 2012, 2016; Eichhorst 2016), demonstrate that there are three species of the genus Neritina registered for the Brazilian Province to date. We emphasize the importance of the separation of N. virginea and N. meleagris in future studies, especially those conducted in the Brazilian Province, since these species may cohabit. In the field, these species can be identified with a detailed observation of the leading edge patterns of their shells, assisting ecological studies. Further research is needed in other areas along the Brazilian Province to determine the geographic distribution of N. virginea and N. meleagris, highlighting the locations where they co-occur.

Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) provided a postdoctoral fellowship to C.X. Barroso (PNPD process number 88882.306440/2018-01) and J.E.P Freitas (PNPD process number 88887.466769/2019-00). H. Matthews-Cascon and T.M.C. Lotufo are research fellows from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). This work was partially funded by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), projects # 2015/17177-6 and 2017/11948-6. The authors would like to thank the Central Analítica-UFC/CT-INFRA/MCTI-SISANO/Pró-Equipamentos CAPES for the support. The authors are also indebted to Patricia N. Bordallo (Embrapa), Vicente V. Faria (UFC), and Diego V. Wilke (UFC) for making laboratory infrastructure available in part of the research. Finally, we are indebted to Dr. Thierry Backeljau, Dr. Ahmed Abdou and the anonymous referees for their careful reviews.