INTRODUCTION
Sympatric speciation remains a controversial topic in evolutionary biology. Several mechanisms could theoretically produce reproductive isolation of sympatric populations resulting in speciation (Friesen et al., Reference Friesen, Smith, Gómez-Díaz, Bolton, Furness, González-Solís and Monteiro2007). These include polyploidization (flowering plants), non-random mating driven by ecological isolation (host preference, differential migratory route, allochrony), and disruptive selection (frequency-dependent competition, divergent sexual selection) (Friesen et al., Reference Friesen, Smith, Gómez-Díaz, Bolton, Furness, González-Solís and Monteiro2007). In marine systems, where swimming larval dispersal is rampant, it is more difficult to envision mechanisms reducing gene flow between sympatric populations. However, a few examples have been identified in marine invertebrates and fish, with causal mechanism including sexual selection (Jones et al., Reference Jones, Moore, Kvarnemo, Walker and Avise2003), feeding specialization (Johannesson et al., Reference Johannesson, Rolán-Álvarez and Ekendahl1995), habitat preference (Stanhope et al., Reference Stanhope, Hartwick and Baillie1993) and host shifts (Munday et al., Reference Munday, van Herwerden and Dudgeon2004).
In this paper we study two putative species of Dondice Marcus, 1958, a poorly-known group of cnidarian-feeding aeolid opisthobranch molluscs with swimming larvae. Because these two species have sympatric ranges but different feeding specializations, this case could constitute an example of sympatric speciation in organisms with larval dispersal.
Dondice is a small group of opisthobranch sea slugs with only three species currently recognized as valid, including the Mediterranean and Eastern Atlantic Dondice banyulensis Portmann & Sandmeier, 1960; the Caribbean and Brazilian species Dondice occidentalis (Engel, 1825); and the highly specialized Dondice parguerensis Brandon and Cutress, Reference Brandon and Cutress1985, which feeds and lives on up-side-down jellies of the genus Cassiopeia. See Edmunds (Reference Edmunds1964) and García & García (Reference García and García1984) for reviews of the systematics of this group.
Dondice parguerensis was originally described from Puerto Rico (La Parguera) as a distinct species from Dondice occidentalis, based on several colour and anatomical differences as well as its very distinct life history, including diet and developmental mode (Brandon & Cutress, Reference Brandon and Cutress1985). Both D. parguerensis and D. occidentalis exhibit a wide range of colour morphs (Figure 1), some of which are very similar between the two species (see Valdés et al., Reference Valdés, Hamann, Behrens and DuPont2006). This makes distinguishing between the two species using external morphology and colour exceedingly difficult. Dondice parguerensis has been reported from Puerto Rico and Panama (Valdés et al., Reference Valdés, Hamann, Behrens and DuPont2006); a range that overlaps with that of D. occidentalis, which is widespread in the Caribbean and Brazil (Valdés et al., Reference Valdés, Hamann, Behrens and DuPont2006). Although not explicitly discussed by Brandon & Cutress (Reference Brandon and Cutress1985), considering D. parguerensis a distinct species from the sympatric D. occidentalis suggests that the feeding specialization of D. parguerensis might have been a predominant behavioural trait acting as a mechanism of genetic isolation and speciation.
Colour variation in western Atlantic Dondice: (A) Dondice parguerensis from Panama (LACM 2004-10.2); (B) Dondice occidentalis from Florida (LACM 2003-41.5); (C–D) Dondice occidentalis from Bahamas.
Freshly collected specimens of D. occidentalis from Florida, Bahamas, and Costa Rica, and D. parguerensis from Puerto Rico (La Parguera) and Panama allowed for the first time to test the validity of these two species using molecular markers. The main objective of this study is to investigate possible molecular differences between the two putative species.
MATERIALS AND METHODS
Source of specimens
Twenty-three specimens of Dondice were sequenced or otherwise examined (Table 1). Specimens were collected by the authors, donated by colleagues, or obtained from museum collections. Specimens collected by the authors were photographed alive, narcotized using a 1 M solution of MgCl2 and preserved in 70% or 100% EtOH. Most specimens are deposited at the Natural History Museum of Los Angeles County, California (LACM). Some small specimens were completely used up for DNA extraction and no tissue samples remain.
List of specimens used in the study; including locality, voucher number, and GenBank Accession numbers.
*, Partial 16S sequence; **, specimens studied morphologically.
Specimens and DNA extraction
DNA extraction was performed using a hot Chelex protocol, a DNeasy Blood and Tissue Kit (Qiagen), or a standardized phenol chloroform extraction technique. Approximately 1–3 mg of tissue was taken from the animals and cut into fine pieces for extraction in each of the three protocols. For the Chelex extraction, the tissue was rinsed and rehydrated using 1.0 ml TE buffer (10 mM Tris, 1 mM EDTA, pH 8.0) for 20 min. A 10% (w/v) Chelex 100 (100–200 mesh, sodium form, Bio-Rad) was prepared using TE buffer. After rehydration, the mixture was then centrifuged, 975.00 ml of the supernatant was removed, and 175.00 ml of the Chelex solution was added. Samples were then heated in a 56°C water bath for 20 min, heated in a 100°C heating block for 8 min, and the supernatant was used for polymerase chain reaction (PCR). The DNeasy protocol supplied by the manufacturer was followed, with some modifications. The elution step was modified such that the first elution was collected using 100.00 ml of Buffer AE and was allowed to incubate at room temperature for 5 min. In a new test tube, a second elution step was conducted using 200.00 ml of Buffer AE and was also allowed to incubate at room temperature for 5 min. The first elution was used for PCR. The phenol chloroform extractions were done by first placing 1–3mg of finely cut tissue into 200 ml of lysis buffer (100 mM EDTA, 10 mM Tris, and 1% SDS), 20 ml of Proteinase K was then added and the tissue samples were allowed to lyse overnight in a heated 56°C water bath. After the enzyme had completely dissolved the tissue, 220 ml of 25:24:1 phenol/chloroform/isoamyl alcohol was added and the solution was vortexed vigorously for 10 s before being centrifuged at maximum speed for 30 s. The top aqueous phase was removed and transferred to a new tube and a volume equalling 1/10 of the supernatant of 3M sodium acetate (pH 5.2) was added and mixed by simple inversion. To this mixture, ice cold 100% ethanol was added in an amount proportional to twice the volume of the solution. The tubes were vortexed and placed in a −80°C freezer for 30 min. Subsequently, the reaction mixtures were centrifuged for 5 min at maximum speed. The supernatant was discarded and 1 ml of 70% ethanol was added, inverted, and centrifuged for another 5 min at maximum speed. Most of the supernatant was then discarded and the remaining pellet was allowed to dry overnight. After drying was complete the pellet was dissolved in 200 ml of TE buffer (pH 8.0) and stored in a −20°C freezer until needed for PCR. The reason why such a variety of extraction protocols were implemented was due to the fact that some tissue samples could not be amplified when using a particular extraction method. It was not uncommon for a particular specimen to be unsuccessfully extracted in all but one of the three methods. The use of the phenol chloroform extraction technique as well as the DNeasy extraction protocol were only used as alternatives to the Chelex extraction when it became apparent that the DNA could not be successfully extracted using this latter technique.
Primers
Colgan's universal H3 primers (Colgan et al., Reference Colgan, McLauchlan, Wilson, Livingston, Edgecombe, Macaranas, Cassis and Gray1998) and Palumbi's universal 16S primers (Palumbi, Reference Palumbi, Hillis, Moritz and Mable1996) were used to amplify the regions of interest for all specimens (Table 2). Additionally, internal primers for 16S designed for another group of opisthobranchs (Ornelas-Gatdula et al., Reference Ornelas-Gatdula, DuPont and Valdés2011) were used, resulting in partial sequences for some specimens. Multiple attempts to obtain CO1 (cytochrome oxidase) sequences using different primers were unsuccessful.
Forward (F) and reverse (R) PCR primers used to amplify regions of the nuclear H3 gene and mitochondrial 16S genes.
PCR amplification and sequencing
The master mix was prepared using 34.75 ul H2O, 5.00 ul Buffer B (ExACTGene, Fisher Scientific), 5.00 ul 25 mM MgCl2, 1.00 ul 40 mM dNTPs, 1.00 ul each 10 mM primer, 0.25 ul 5 U/ul Taq, and 2.00 ul extracted DNA. Reaction conditions for H3 and 16S rRNA were as follows: an initial denaturation for 2 min at 94°C, 35 cycles of (1) denaturation for 30 s at 94°C, (2) annealing for 30 s at 50°C, and (3) elongation for 1 min at 72°C, and a final elongation for 7 min at 72°C. PCR products yielding bands of appropriate size were purified using the Montage PCR Cleanup Kit (Millipore). Cleaned PCR samples were quantified using a NanoDrop 1000 Spectrophotometer (Thermo Scientific). Samples were sequenced at the City of Hope DNA Sequencing Laboratory (Duarte, California) using chemistry type BigDye V1.1.
Phylogenetics
Sequences for each gene were assembled and edited using Geneious Pro 4.7.4 (Biomatters Ltd.). Geneious was also used to extract the consensus sequence and to construct the alignment for each gene using the default parameters. The sequences were not trimmed after alignment. A total of 424 bp for 16S (232 bp in the partial sequences) were used for the phylogenetic analyses. An analysis was conducted for 16S only because H3 was not phylogenetically informative (only one allele was present in the ingroup). For outgroup comparison we selected the species Godiva quadricolor (Barnard, 1927) using sequences available in GenBank (HM162680). Dondice has been considered a synonym of Godiva and these two groups are morphologically similar (Edmunds, Reference Edmunds1964; García & García, Reference García and García1984) suggesting a close evolutionary history. Sequences of Dondice banyulensis Portmann & Sandmeier, 1960 from the Mediterranean were obtained from GenBank (HQ616740, GQ403751) and included in the analysis for comparison. A second analysis was conducted excluding the sequence HQ616740 due to misidentification of that specimen (see Discussion).
Maximum likelihood analyses were conducted using PAUP*4.0 (Swofford, Reference Swofford2002). Robustness of each clade was assessed by bootstrap support (Felsenstein, Reference Felsenstein1985) based on 2000 replicates with heuristic search, TBR branch-swapping algorithm, and 100 random additions.
Bayesian analyses were executed in MrBayes v3.1.2 (Huelsenbeck & Ronquist, Reference Huelsenbeck and Ronquist2001). The Markov chain Monte Carlo analysis was run with two runs of six chains for five million generations, with sampling every 1000 generations. Effective sample sizes and convergence of runs were assessed using Tracer 1.4.5 (Rambaut & Drummond, Reference Rambaut and Drummond2009). The default 25% burn-in was applied before constructing majority-rule consensus tree.
The Akaike information criterion (Akaike, Reference Akaike1974) was executed in jModelTest (Posada, Reference Posada2008) and MrModeltest v2.3 (Nylander, Reference Nylander2004), to determine the best-fit model of evolution.
Population genetics
The population structure and the genetic differentiation between populations of western Atlantic Dondice were analysed for the 16S gene, using analysis of molecular variance (AMOVA) and the FST analysis, both implemented in Arlequin 3.5 (Excoffier & Lischer, Reference Excoffier and Lischer2010). Significance of the AMOVA and FST analyses was tested using 16,000 permutations. Samples were assigned into two groups according to their diet: ‘Parguerensis’ (Dondice parguerensis including animals from Puerto Rico and Panama) and ‘Occidentalis’; (D. occidentalis). For the AMOVA, the proportion of variation within populations (ΦST) and between populations (ΦSC) were computed.
Radulae
In order to test whether feeding preference had influenced the evolution of radular morphology in Dondice, the buccal masses of several specimens (Table 1) were dissected under a light microscope (Nikon SMZ-100). Each buccal mass was dissolved in 10% sodium hydroxide until the radula and jaw were isolated from the surrounding tissue. The radula and jaw were then rinsed in water, dried, mounted and sputter coated for examination with a scanning electron microscope (SEM) Hitachi S-3000N at the Natural History Museum of Los Angeles County.
RESULTS
Molecular data
All H3 sequences of western Atlantic Dondice were identical, representing the same allele. This includes specimens of D. parguerensis and D. occidentalis. For the 16S gene, five different haplotypes were identified in western Atlantic Dondice, uncorrected pairwise identity = 99.6%. All specimens of D. parguerensis had two unique substitutions not found in D. occidentalis.
The analysis of the 16S gene produced similar Bayesian and maximum likelihood consensus trees. In both trees all western Atlantic specimens of Dondice are clustered into a well-supported monophyletic group (posterior probability = 1, bootstrap support = 99). Within this group there is a poorly supported second clade containing all specimens identified as D. parguerensis (posterior probability = 0.82, bootstrap support = 61). There is no structure or support for the monophyly of specimens identified as D. occidentalis. A second set of analyses excluding the problematic sequence HQ616740 (see Discussion) provided similar trees but with higher support values (Figure 2). In Bayesian and maximum likelihood consensus trees all western Atlantic specimens of Dondice are clustered into a well-supported monophyletic group (posterior probability = 1, bootstrap support = 100). Within this group there is a second clade containing specimens identified as D. parguerensis (posterior probability = 0.94, bootstrap support = 68). There is no structure or support for the monophyly specimens identified as D. occidentalis except for a poorly supported clade containing some specimens from the Bahamas (posterior probability = 0.65).
Bayesian consensus tree for 16S. Support values for likelihood and Bayesian approaches are shown only for nodes with bootstrap values greater than 50 or posterior probabilities greater than 0.5.
The FST analysis of 16S haplotypes revealed that D. occidentalis and D. parguerensis are significantly distinct (FST = 0.40669, P = 0.00548 ± 0.0018). For the AMOVA, among population differentiation explained 40.67% of the covariance component, while within-population differentiation explained 59.33% of the covariance component (ΦST = 0.40669, P = 0.00631 ± 0.0006).
Morphological data
The radular formula of the specimens examined varies from 13 × 0.1.0 to 22 × 0.1.0. All specimens have rachidian teeth with a central cusp surrounded by 5–7 sharp denticles on each side (Figure 3). In some specimens the innermost denticle emerges from the base on the central cusp rather than the base of the teeth. The rachidian teeth are similar in shape and size throughout the radula. Comparison of the radulae of both groups of western Atlantic Dondice shows no consistent pattern. There is no correlation between tooth morphology and diet. For example, specimens found on Cassiopeia or on hydroids can have denticles emerging from the base of the central cusp or not.
Scanning electron micrographs of the radulae of Caribbean specimens of Dondice: (A–B) Dondice parguerensis, Panama (LACM2004-10.2): (A) anterior teeth and (B) posterior teeth; (C–D) Dondice occidentalis, Costa Rica (INB0003348310): (C) anterior teeth and (D) posterior teeth; (E–F) Dondice occidentalis, Bahamas (LACM 177708): (E) anterior teeth and (F) posterior teeth; (G–H) Dondice occidentalis, Bahamas (LACM 177709): (G) anterior teeth and (H) posterior teeth; (I–J) Dondice occidentalis, Florida (LACM 2003-41.5): (I) anterior teeth and (J) posterior teeth.
DISCUSSION
The results of the analyses conducted in this study are inconclusive as to whether Dondice occidentalis and D. parguerensis are genetically distinct. The two putative species are identical in the H3 gene, which is generally more conserved and in some cases too conserved to recover species level differences (Ornelas-Gatdula et al., Reference Ornelas-Gatdula, DuPont and Valdés2011). A H3 sequence of D. banyulensis available in GenBank (HQ616804) has 34 substitutions when compared with D. occidentalis/parguerensis, which would indicate that H3 provides ample resolution to distinguish species within Dondice. However, this GenBank sequence (HQ616804) is potentially problematic. It was submitted along with a 16S sequence of the same specimen (HQ616740), also included in the present analysis, which is substantially different (90 substitutions) from another 16S sequence of D. banyulensis also available in GenBank (GQ403751). BLASTing the problematic sequence (HQ616740) reveals that it is identical to a sequence of the aeolid nudibranch Dicata odhneri Schmekel, 1967. Thus we suspect this sequence was obtained from a misidentified animal or it was mistakenly uploaded by the authors in GenBank. These two problematic sequences were not included in the final analysis presented here (Figure 2).
Contrary to the H3 results, all specimens of D. parguerensis from two geographically distant locations (Puerto Rico and Panama) have unique substitutions in the 16S gene (genetic synapomorphies) that suggest reproductive isolation from D. occidentalis. The 16S gene is variable enough to distinguish among species of opisthobranchs of diverse clades (Turner & Wilson, Reference Turner and Wilson2007; Anthes et al., Reference Anthes, Schulenburg and Michiels2008; Johnson & Gosliner, Reference Johnson and Gosliner2012, Krug et al., Reference Krug, Händeler and Vendetti2012a) and is even informative in population genetics studies (Krug et al., Reference Krug, Asif, Baeza, Morley, Blom and Gosliner2012b). Moreover, the FST analysis of 16S confirmed that the two species (D. parguerensis and D. occidentalis) are genetically distinct. Conversely, the Bayesian and maximum likelihood analyses do not recover D. occidentalis and D. parguerensis as reciprocally monophyletic, and the AMOVA results are inconclusive, questioning the separation of the two putative species.
Although we were unable to find any consistent radular differences between D. occidentalis and D. parguerensis, Brandon & Cutress (Reference Brandon and Cutress1985) described several differences in coloration, size, number of ceratal groups and more importantly larval development. According to these authors D. occidentalis has lecitotrophic development whereas D. parguerensis has planktotrophic development. We have observed differences in coloration, with specimens of D. parguerensis never having any red markings. Because of these morphological differences and the presence of unique substitutions in 16S gene, we prefer to maintain D. occidentalis and D. parguerensis as different species.
One possible explanation for the inconclusiveness of the molecular results is that D. occidentalis and D. parguerensis are in the process of diverging into two different species but the speciation process is incipient or incomplete. Examination of additional genes and/or microsatellite data as well as mating preference assays could provide additional insight into assortative mating patterns within or between the two populations/species helping to resolve the question of how many species are there. This is beyond the scope of this paper.
Regardless of the completeness of the speciation process, the data presented here suggest that feeding specialization of some Dondice populations on a novel food source (up-side-down jellies) has possibly lead to the reduction of gene flow with hydroid-eating populations. The fact that distant populations of D. parguerensis found in Panama and Puerto Rico have the same substitutions in the 16S gene, not present in any D. occidentalis, suggests that genetic isolation is maintained across the range of the two putative species and represents a potential example of sympatric speciation in a benthic marine organism. Although sympatric speciation has been proposed for other aquatic and marine organisms, it is often difficult to provide conclusive evidence or completely reject alternative hypotheses (e.g. Barluenga et al., Reference Barluenga, Stölting, Salzburger, Muschick and Meyer2006; Schliewen et al., Reference Schliewen, Kocher, McKaye, Seehausen and Tautz2006). Bolnick & Fitzpatrick (Reference Bolnick and Fitzpatrick2007) commented that in oceanic habitats where range expansions are easier, case studies of sympatric speciation might rest on ecological evidence, rather than on biogeographic evidence against geographic isolation. Host shifts are expected to occur within the range of a population, and can directly confer reproductive isolation if mating occurs on the host (Bolnick & Fitzpatrick, Reference Bolnick and Fitzpatrick2007). Consequently, sympatric sister species using different hosts provide a reliable line of evidence for sympatric speciation. For example, strong host fidelity along with disruptive selection has been proposed as a mechanism of rapid reproductive isolation resulting in sympatric speciation in Gobiodon fish (Munday et al., Reference Munday, van Herwerden and Dudgeon2004). In the case of Dondice parguerensis adult individuals and egg masses are found exclusively in up-side-down jellies (Brandon & Cutress, Reference Brandon and Cutress1985; A.V., personal observation), which is most likely a host shift from the ancestral hydroid feeding condition and a possible mechanism of reproductive isolation from D. occidentalis. Other examples of dietary specialization acting along with disruptive selection have been proposed as mechanisms of reproductive isolation resulting in sympatric speciation in marine molluscs and fish (Johannesson et al., Reference Johannesson, Rolán-Álvarez and Ekendahl1995; Munday et al., Reference Munday, van Herwerden and Dudgeon2004). Further genetic and behavioural research would be valuable in determining if and/or how ecological host-shifting has driven sympatric speciation in these two Dondice species.
ACKNOWLEDGEMENTS
We are extremely grateful to Anne DuPont and Patrick Krug for providing specimens essential for this project. Lindsey Groves from the LACM helped with access to the collection and curation of specimens. Yolanda Camacho kindly made available specimens from the Universidad de Costa Rica. An anonymous referee made constructive comments that greatly improved this manuscript.
FINANCIAL SUPPORT
This paper was partially supported by the National Science Foundation-funded PEET project ‘Phylogenetic Systematics of Nudibranchia’ (DEB-0329054). The SEM work was conducted at the Natural History Museum of Los Angeles County, supported by the National Science Foundation under MRI grant DBI-0216506 with the assistance of Giar-Ann Kung. Fieldwork in Panama was supported by a postdoctoral scholarship of the Smithsonian Institution Tropical Research Institute to the junior author.
