Abstract
Molecular and morphological evidence indicates that the widespread North Pacific species Diaulula sandiegensis is a complex of two pseudocryptic species. One ranges from the Sea of Japan to northern California, whereas the other is found from northern California to Baja California. The name Diaulula sandiegensis applies to the southern species and Diaulula odonoghuei Steinberg, 1963 to the northern species; a neotype is designated for the latter. The two have consistent differences in colour pattern that were overlooked in previous studies. Field observations from San Diego, California, to Bamfield, British Columbia, revealed that the ranges of the two species partially overlap: D. odonoghuei is restricted to intertidal and bay habitats from Fort Bragg, California northwards, whereas D. sandiegensis is found throughout the entire region and is more common in subtidal habitats. Both species feed primarily on a sponge, Haliclona sp., in intertidal habitats, but faecal analysis indicated that D. sandiegensis in subtidal habitats feeds on Neopetrosia problematica. Haplotype-network analyses reveal substantial levels of genetic diversity in both species, which is highly structured in D. odonoghuei, suggesting a complex history for this species. Field surveys, a common garden experiment and a mating study were conducted to investigate the presence of assortative mating. During two mating trials, individuals that encountered one another for the first time always mated with the same species (11 matings). Reproductive isolation in the form of premating barriers may explain the absence of gene flow between these geographically overlapping sister species. A molecular-clock analysis indicates that glacial events may have played a role in speciation and the origin of population structure in the D. sandiegensis species complex, possibly through climate-driven vicariant events.
INTRODUCTION
The North Pacific Ocean is a region of great interest for biogeographers and marine evolutionary biologists. Pleistocene glacial cycles have resulted in repeated range expansions and contractions in temperate and boreal species in the Northern Hemisphere (Hewitt, 2000). The formation of ice sheets in large areas of the Northern Hemisphere during glacial periods and significant climate cooling resulted in the local extinction of numerous species as coastal areas became uninhabitable due to accumulation of ice along the shoreline (Dyke et al., 2002). When the ice sheets retreated during interglacial periods, populations recovered with influx of larvae or adults from refugia in southern latitudes, allowing a rapid northward expansion during the transition from the Pleistocene to the Holocene (Marko, 1998; Hellberg, Balch & Roy, 2001; Edmands & Harrison, 2003). This process has been documented in eastern Pacific marine benthic organisms (Jacobs, Haney & Louie, 2004), including heterobranch sea slugs (Ellingson & Krug, 2006; Cooke et al., 2014).
A second possible effect of glacial cycles is the generation of biodiversity. Speciation can occur allopatrically following vicariance between populations isolated in refugia by ice formation or changes in sea level (Liu et al., 2011; Shen et al., 2011), or sympatrically as a result of adaptive radiation during range expansion in interglacial periods (Schluter & Rambaut, 1996; Pigeon, Chouinard & Bernatchez, 1997). Sibling-species pairs arising from speciation driven by Pleistocene glaciations should display low levels of genetic divergence.
In the North Pacific Ocean, several groups of heterobranch sea slugs include closely related but distinct species that occur in the western and eastern Pacific; however, some species have unusually large geographic ranges across the region. For example, Diaulula sandiegensis (Cooper, 1863) is a widespread nudibranch recorded from Mexico to Japan (Behrens & Hermosillo, 2005) and Korea (Koh, 2006), with broad variation in colour pattern (Behrens & Valdés, 2001). Several studies attempted to explain the biological nature of this colour variation. McDonald (1983) indicated that specimens of D. sandiegensis from bays are darker than those from open coastal waters. Penney (2013a) recognized at least two sympatric colour forms of D. sandiegensis in British Columbia, which differed in the shape of the dark markings. He identified a ‘ring’ form with few brown rings symmetrically distributed and offset from the midline, typically on a light background, and a ‘spotted’ form with multiple, irregular, dark patches similar to leopard spots in that they are surrounded by a ring of different (in this case, lighter) colour, which often grade into this colour towards the centre. Penney (2013a) also found consistent biological differences between the two forms, for example ‘spotted’ specimens were significantly smaller than ‘ring’ specimens, and whereas ‘spotted’ specimens were predominantly intertidal, ‘ring’ specimens were exclusively subtidal. He argued that these differences in colour, habitat and size could indicate the existence of two distinct species. In contrast, Behrens & Valdés (2001) studied the anatomy of D. sandiegensis across most of its range in North America and recognized a single species. More recently, Kelly (2013) investigated the existence of two species within D. sandiegensis in Humboldt Bay, California, by examining morphological differences, bathymetric range, breeding behaviour and feeding preference. She identified two morphs, the ‘leopard’ morph, with spotting extending onto the mantle margin and the ‘ring’ morph without spots on the margin. The reported diet of D. sandiegensis consists of four to five different species of sponges (Elvin, 1976; Penney, 2013b), but Kelly (2013) indicated that in Humboldt Bay the ‘ring’ morph feeds on Neopetrosia problematica, a yellow to white encrusting sponge, whereas the ‘leopard’ morph feeds on Haliclona permollis, a purple encrusting sponge. She speculated that this difference in dietary preference was evidence to support the existence of two distinct species within D. sandiegensis.
To further address whether there are one or more distinct species within nominal D. sandiegensis in the North Pacific Ocean, we studied its genetic structure across most of its range including populations from the northeastern Pacific (Russia and Japan). We also studied mating behaviour, habitat use and food preferences of the two colour forms in areas of the northwestern Pacific where they co-occur. Finally, we quantified variation in spotting pattern and radular morphology in the two forms.
MATERIAL AND METHODS
Source of specimens
Diaulula sandiegensis specimens were obtained by SCUBA diving, on floating docks, from the shoreline during low tide, or were donated by colleagues. A specimen of D. hummelincki (Ev. Marcus & Er. Marcus, 1963) from the Caribbean Sea, used as the outgroup in molecular phylogenetic analyses, was obtained from the California State Polytechnic University Invertebrate Collection (CPIC). Specimens of D. greeleyi (MacFarland, 1909) and D. nayarita (Ortea & Llera, 1981), collected in the Caribbean and eastern Pacific respectively, were used for calibration of molecular-clock analyses; these were obtained from the Natural History Museum of Los Angeles County (LACM). Specimens collected by the authors were photographed, preserved in 95% ethanol and deposited in the CPIC.
Morphological variation and phenotypic plasticity in spotting pattern
A total of 444 specimens of D. sandiegensis were photographed in situ at 14 field sites from San Diego, California, to Bamfield, British Columbia, in 2009, 2010 and 2011. Additionally, 88 D. sandiegensis images were obtained, with permission, from the World Wide Web. These images were analysed for variation in dorsal spotting pattern and dorsal ‘background’ coloration, in order to distinguish two different morphotypes: the spotted morph and the ringed morph (corresponding to the ‘leopard’ and ‘ring’ morphs, respectively, of Kelly, 2013). Each specimen was photographed with a ruler for scale. ImageJ software (Abramoff, Magalhães & Ram, 2004) was used to measure the rhinophore-to-gill length (RGL; from the centre of the rhinophores to the centre of the gill) as well as the number and diameter of all spots on the dorsum. RGL was used to represent total length, so that nudibranchs could be measured from photographs without removing them from their natural substrate. This technique also avoids problems due to curling of the mantle margin, which can obscure the true length. Photographs of 321 specimens of D. sandiegensis were analysed in this way, to distinguish the morphotypes and to determine the overlap in colour pattern between them (in terms of number, type and arrangement of dorsal spots as well as background coloration of the dorsum). RGL was obtained from the photographs of 233 individuals. Since Behrens & Valdés (2001) found a wide variation in the number of dorsal spots on living D. sandiegensis from San Diego, California, to Vancouver Island, British Columbia, a linear regression was used to determine whether the number of dorsal spots was dependent on RGL. A t-test on log-transformed data was used to determine whether individuals with and without spots on the mantle margin differed in the number of dorsal spots.
Radulae of ten specimens were prepared as follows. Forceps were used to open the dorsum between the rhinophores and the buccal mass extracted through this incision, then placed in 10% sodium hydroxide solution for 3 d to dissolve the tissue. Radulae were mounted and images taken at the LACM scanning electron microscopy laboratory using a Hitachi S-3000N variable-pressure scanning electron microscope (SEM).
To determine whether spotting pattern is plastic, 36 specimens (32 spotted and 4 ringed morphs) were located and photographed repeatedly in the field at three sites at Trinidad, California, over the course of 2–10 weeks. Individual animals were identified in the field by their unique dorsal spotting patterns. Individuals were tracked daily over 3–4 d during a series of low tides and then located again in the same general area, c. 2 weeks later, during the next low-tide cycle. For each photograph, RGL and spot sizes were measured and all dorsal spots were counted using ImageJ. To further investigate changes in dorsal spotting pattern with growth, ‘common garden’ experiments were conducted at the Telonicher Marine Laboratory (Trinidad, California). Twenty-eight individuals were collected in 2010: 14 from the subtidal of Monterey Bay at Shale Reef (SCUBA diving at 15 m depth) and 14 from the intertidal of Endert‘s Beach at Crescent City, California (intertidal collection). All individuals found in the Crescent City intertidal had the pattern of the spotted morph and all individuals found in the Monterey Bay subtidal had the pattern of the ringed morph. Each nudibranch was kept separately for 24 h in a plastic container with a fine-mesh screen lid to allow collection of faeces. The spicules of sponges consumed in the field were identified in the faeces using light and SEM images following Lee, Elvin & Reiswig (2007). Animals were kept in an aquarium with a closed seawater system; seawater was aerated and kept at 11 °C. The two field-collected populations were fed the same food regime, an ad libitum supply of the sponge Haliclona ‘sp. A’ (Lee et al., 2007), collected from floating docks in Humboldt Bay, California. Individuals were repeatedly photographed over the course of 18 weeks and images analysed as described above.
Latitudinal distribution and habitat use
The images from the 444 specimens of D. sandiegensis found during field investigations in 2009, 2010 and 2011 and the 88 images from the World Wide Web (see above) provided a photographic record of morphotypic frequencies across 25 sites from San Diego, California, to Barkley Sound, British Columbia.
Variation in morphotypes across different habitat types was observed by Morris et al. (1980), McDonald (1983) and Penney (2013a). Therefore, morphotypic frequencies were compared across the same three habitat types used by these authors: (1) coastal subtidal sites, (2) coastal intertidal habitats (mostly rocky shores) and (3) bay habitats.
Sponge prey
To identify sponge prey associated with D. sandiegensis in the intertidal zone, a field investigation of nudibranch–sponge associations was conducted in 2009 using 178 animals found in intertidal habitats from San Diego, California to Barkley Sound, British Columbia. Specimens were identified as feeding if they were found in proximity to sponges that had sections or paths removed, the width of which was similar to the width of the nudibranch (see Supplementary Material Fig. S1A for an example). Since individual nudibranchs could be identified by their unique dorsal spotting pattern, it was possible to watch the progress of sponge removal over several days by the same nudibranch. A Pearson Chi-square test was used to compare sponge use by the two morphs. Sponge spicules and skeletal structure (from pieces of sponge obtained in the field) were examined with light and SEM to confirm field identifications (Lee et al., 2007). Spicules from prey sponges remain intact during digestion and can be found in the faeces of dorid nudibranchs (Bloom, 1981). To identify sponge prey used by D. sandiegensis in a subtidal habitat, a faecal analysis was conducted using the 14 D. sandiegensis collected from Monterey Bay (see above).
Reproductive compatibility
To investigate the reproductive compatibility of the morphs in their natural environment, mating events were documented from 44 photographs of mating pairs taken in the field in 2009, 2010 and 2011. Like other nudibranchs, D. sandiegensis is a simultaneous hermaphrodite and mates reciprocally. We defined two individuals as mating when they were positioned right side to right side, stationary and with the mantle margin of each pushed up to expose the white underside of the mantle. For these mating events, captured in photographs, it was impossible to confirm sperm transfer. However, animals in this position always showed copulation with at least one penis uniting them when they were pulled apart (Kelly, 2013).
In addition, two mate-choice trials were conducted at the Telonicher Marine Laboratory to determine whether like morphs would mate preferentially. The first of these trials included ten specimens collected in October 2010, from three sites: (1) Crescent City, California, intertidal zone; (2) Trinidad, California, intertidal zone and (3) floating docks in Humboldt Bay. These ten specimens (six spotted and four ringed morphs) were placed together in one aquarium with an open circulating seawater system, allowing each animal to choose mates at random. The initial contact, mating events, duration of each mating event, morphotype of each partner and premating behaviour were documented using sequences of photographs captured every 60 s with a Nikon Coolpix® 995 camera attached to a Digisnap 2000® intervalometer (Harbortronics). In this study, a mating event between unique individuals was defined as the first photographed mating for at least one individual in the pair. Individual animals were identified by their unique dorsal spotting patterns. After documenting these initial mate choices, recording was continued in three subsequent runs (total 87.5 h) to determine whether these individuals mated repeatedly. The first run lasted 14.3 h; the second, using nine of the same specimens (minus one spotted morph that died), lasted 13.8 h; the third, using six from the original group (minus two spotted and one ringed morphs that died) lasted 59.4 h. If food was provided during a trial, the nudibranchs would spend most of their time feeding rather than mating; therefore, between runs each was held in a separate container and fed an ad libitum supply of Haliclona sp. A for 1–2 weeks. The run times differed due to adjustments of pixel quality on the camera used to allow more recording time.
The second mating trial used eight new D. sandiegensis (four spotted and four ringed morphs) collected from the same sites as the first, in July 2011. Similarly, after the initial mate choices were documented, behaviour was recorded for a total of 232 h in four runs. The first run lasted 71 h; the second used the same animals plus one additional spotted morph from Humboldt Bay and ran for 42.4 h; the third used eight of these same animals (minus one spotted morph that died) and lasted 118.6 h. A fourth trial, using these same eight animals, was conducted in August 2011, to observe mating among specific morph combinations. Four pairs of dissimilar morphs were placed in individual aquaria and observed for 22.3 h.
DNA extraction, amplification and sequencing
Map showing the localities sampled for molecular analysis of the Diaulula sandiegensis species complex. The circle sizes are proportional to the number of specimens sequenced; black sector indicates ringed morph (D. sandiegensis s. s.); white indicates spotted morph (D. odonoghuei).
DNA extraction was performed using DNeasy Blood and Tissue Kit (Qiagen) or a standard Chelex® extraction. A 1 mm3 piece of foot or mantle tissue was used. For Chelex extraction, the macerated tissue was transferred into a 1.7-ml microcentrifuge tube containing 1 ml of 1× TE buffer and agitated for 20 min to rehydrate the tissue and begin to disassociate the cells. Samples were then vortexed for 5 s, centrifuged at 23,897 × g for 3 min and 975 μl of supernatant removed from each. 175 μl of 10% Chelex solution was added to each sample, vortexed and samples placed in a water bath at 56 °C for 20 min. After vortexing for 5 s they were placed in a heat block at 100 °C for 8 min, vortexed again for 5 s and then centrifuged at 23,897 × g for 3 min. The resulting supernatant was used for PCR. For DNeasy extraction, the manufacturer‘s protocol for tissue samples was followed.
The polymerase chain reaction (PCR) was used to amplify portions of the mitochondrial cytochrome c oxidase I (COI) and 16S ribosomal RNA (16S) genes, as well as the nuclear histone H3 (H3) gene and the first 500 bp of the 18S ribosomal RNA (18S) gene. The following universal primers were used: COI (LCOI490 5′-GGTCAACAAATCATAAAGATATTGG-3′, HCO2198 5′-TAAACTTCAGGGTGACCAAAAAATCA-3′; Folmer et al., 1994), 16S rRNA (16Sar-L 5′-CGCCTGTTTATCAAAAACAT-3′, 16Sbr-H 5′-CCGGTCTGAACTCAGATCACGT-3′; Palumbi, 1996), H3 (H3 AF 5′-ATGGCTCGTACCAAGCAGACGGC-3′, H3 AR 5′-ATATCCTTGGGCATGATGGTGAC-3′; Colgan et al., 1998) and 18S partial (18SA1 5′-CTGGTTGATCCTGCCACTCATATGC-3′, 18S700R 5′-CGCGGCTGCTGGCACCAGAC-3′). Amplification of DNA was confirmed using agarose gel electrophoresis with ethidium bromide. PCR products were sent to Source Bioscience Inc. (Santa Fe Springs, CA, USA) for sequencing. Sequences were assembled, edited and aligned using Geneious Pro v. 8.1 (http://www.geneious.com; Kearse et al., 2012).
Phylogenetic analyses and haplotype network
The best-fit models of evolution for each gene were determined by using the Akaike information criterion (Akaike, 1974) implemented in jModelTest (Posada, 2008); these were HKY+G for COI and HKY+I for 16S; H3 and 18S were not included in the analyses because all sequences were identical and therefore not informative. Maximum-likelihood (ML) analyses were conducted for the two mitochondrial genes (COI, 16S) concatenated with RAxML (Stamatakis, 2006) with 10,000 bootstrap replicates and the GTRGAMMAI model (GTR+I+G). Bayesian analyses were run using MrBayes v. 3.2.1 (Ronquist & Huelsenbeck, 2003), partitioned by gene (unlinked), with two runs of six chains for 10 million generations, sampling every 1,000 generations. Effective sample sizes and convergence of runs were assessed using Tracer v. 1.4.1 (Rambaut & Drummond, 2007) and the default 25% burn-in was applied before constructing majority-rule consensus trees. Bayesian posterior probability (PP) values greater than or equal to 0.95 (Alfaro, Zoller & Lutzoni, 2003) and ML bootstrap (BS) values greater than or equal to 70 (Hillis & Bull, 1993) were considered significant.
A haplotype network was constructed for COI using TCS v. 1.21 (Clement, Posada & Crandall, 2000), with the connection limit set to 95% (Hart & Sunday, 2007). The same program was used to determine relationships between haplotype networks and the geographical location of haplotype groups.
Automatic barcode gap discovery analysis
Automatic barcode gap discovery (ABGD) analysis (Puillandre et al., 2012) was run on the COI ingroup sequences to compare its delimitation of species with the groups identified by the phylogenetic and morphological analyses. ABGD infers the number of species present in a set of sequence data based on gaps in the distribution of pairwise distances between each sequence (Puillandre et al., 2012). The analysis was run for each gene individually using Kimura 2-parameter (K2P) and Tamura–Nei (TN) distance matrices, as well as an uncorrected matrix to prevent bias introduced by use of corrected distances (Srivathsan & Meier, 2011; Collins & Cruickshank, 2012). The matrices were loaded into the online ABGD webtool (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html). The default relative gap width (X) of 1.5 and a range of prior values for maximum divergence of intraspecific diversity (P) from 0.001 to 0.1 were used.
Molecular-clock analysis
A molecular-clock analysis was run on the Bayesian phylogenetic tree based on concatenated mitochondrial sequences of D. sandiegensis. The clock was calibrated using the node between D. greeleyi and D. nayarita set at 3.5 Ma, based on the date of final closure of the Isthmus of Panama (Coates & Obando, 1996) that separated these two likely sister species (Alvim & Pimenta, 2013), and clock estimates for other groups of nudibranchs and sea slugs (Churchill, Valdés & Ó Foighil, 2014; Goodheart et al., 2015). Divergence times were estimated using r8s v. 1.5 (Drummond & Rambaut, 2007), using a combination of one of three methods (Langley–Fitch, non parametric rate smoothing or penalized likelihood) and one of four algorithms (local molecular clock, Powell, quasi-Newton or truncated Newton). Substitution rates per unit of time were obtained for COI using the Langley–Fitch method in order to compare the results with other studies.
RESULTS
Morphological variation and phenotypic plasticity in spotting pattern
Photographs of live animals from across most of the range of Diaulula sandiegensis species complex showing variation in colour pattern, classified into two colour morphs: spotted (A, B, E–G) and ringed (C, D, H). Note lack of spotting on mantle margin in (C,D, H). A, B, E–G.D. odonohuei. C, D, H.D. sandiegensis s. s.A. Russia (CPIC–01264). B. Japan (TL267). C. California (CPIC–00911). D. California (CPIC–00910). E. Japan (CPIC–01263). F. Alaska (CPIC–01072). G. Washington (CPIC–00913). H. California (CPIC–00888). Scale bars = 1 mm.
In the region of the northeastern Pacific where both morphs occur, the field survey of 233 individuals revealed that 99% of the 169 individuals identified as belonging to the spotted morph had solid dorsal spots. In contrast, 76% of the 64 individuals identified as belonging to the ringed morph had ring-shaped spots, while 24% had solid spots (Table 1). This conclusion was guided by the results of the molecular analyses (see below). Thus, the diagnostic trait separating the two morphs is not the shape of the spots, but rather the presence or absence of spots on the mantle margin: all of the ringed morphs (even those with solid spots) can be distinguished from the spotted morph by the absence of spots on the mantle margin. In addition, the number of spots in the spotted morph varied from 23 to 234 (mean 78), whereas in the ringed morph this varied from 4 to 41 (mean 10). The range of overlap (23–41 dorsal spots) included 9% of the spotted morphs and 1% of the ringed morphs. These 9% (15 of 169) of spotted morphs within this zone of overlap were also the smallest individuals, with RGL of less than 1.2 cm (Supplementary Material Fig. S2); they also had clear spots on the mantle margin. Only three ringed morphs fell within this overlap zone and all lacked spots on the mantle margin (Supplementary Material Figs S1C, S2). Overall, individuals with spots on the mantle margin had significantly more spots than those without (t = −25.80, df = 138, P < 0.005, n = 321). Therefore, the presence or absence of spots on the mantle margin distinguished the two morphs over all sizes examined (0.2–8.3 cm).
Description of spotting pattern variability within Diaulula sandiegensis s. l., based on field photographs of 233 animals (169 spotted and 64 ringed morphs). Distinguishing feature of the two morphs is presence or absence of mantle margin spots in spotted and ringed morphs, respectively.
| Mantle margin spots | Spot number | Spot shape | Outer light ring |
|---|---|---|---|
| Present | 23–234 | Solid (99%) | Present (84%) |
| Absent (16%) | |||
| Ring (1%) | Absent | ||
| Absent | 4–41 | Ring (76%) | Absent (99%) |
| Present (1%) | |||
| Solid (24%) | Absent |
| Mantle margin spots | Spot number | Spot shape | Outer light ring |
|---|---|---|---|
| Present | 23–234 | Solid (99%) | Present (84%) |
| Absent (16%) | |||
| Ring (1%) | Absent | ||
| Absent | 4–41 | Ring (76%) | Absent (99%) |
| Present (1%) | |||
| Solid (24%) | Absent |
Description of spotting pattern variability within Diaulula sandiegensis s. l., based on field photographs of 233 animals (169 spotted and 64 ringed morphs). Distinguishing feature of the two morphs is presence or absence of mantle margin spots in spotted and ringed morphs, respectively.
| Mantle margin spots | Spot number | Spot shape | Outer light ring |
|---|---|---|---|
| Present | 23–234 | Solid (99%) | Present (84%) |
| Absent (16%) | |||
| Ring (1%) | Absent | ||
| Absent | 4–41 | Ring (76%) | Absent (99%) |
| Present (1%) | |||
| Solid (24%) | Absent |
| Mantle margin spots | Spot number | Spot shape | Outer light ring |
|---|---|---|---|
| Present | 23–234 | Solid (99%) | Present (84%) |
| Absent (16%) | |||
| Ring (1%) | Absent | ||
| Absent | 4–41 | Ring (76%) | Absent (99%) |
| Present (1%) | |||
| Solid (24%) | Absent |
The number of dorsal spots was correlated with RGL in the spotted morph (R² = 0.48, n = 169, P < 0.0005), but not in the ringed morph (R² = 0.035, n = 64, P = 0.137) (Supplementary Material Fig. S2).
Dorsal background colour varied widely, but the spotted morph tended to have a darker colour than the ringed morph. Of the 175 spotted morph specimens, 50% had a dark brown background colour, 35% light brown, 1% orange and 15% white. Of the 146 ringed morphs, 56% had a white background colour, 30% light brown, 9% orange and 5% dark brown.
Results of the laboratory ‘common garden’ experiment and field survey revealed that all 46 of the spotted morph and all 18 of the ringed morph maintained their spot type (solid or ring-shaped) as they grew, regardless of experimental conditions. Dorsal spots enlarged with growth in both morphs. As spotted morphs grew, additional spots appeared on the edge of the mantle, to maintain spots over the whole dorsal surface. In contrast, the number of spots for all 18 of the ringed morph, from the field as well as the laboratory, remained the same regardless of increasing body length and were never present at the mantle margin.
The radula formula showed no consistent differences between the morphs, varying between 21 × (21–22.0.21–22) for the ringed morph and 20–23 × (29–35.0.29–35) for the spotted morph. Also, no consistent differences were observed between the two groups in tooth morphology (Supplementary Material Figs S3, S4) or the size of the radula (in agreement with Behrens & Valdés, 2001).
Latitudinal distribution and habitat use
Locations of field sites where 444 individuals of Diaulula sandiegensis s. l. were photographed in California, Oregon, Washington and British Columbia, and of locations of 88 photographs obtained from the World Wide Web (asterisks). Numbers in parenthesis on the left are sample sizes for each environment (coastal subtidal, coastal intertidal, and bay). Pie charts indicate frequencies of morphotypes. Ringed morph corresponds to D. sandiegensis s. s., spotted morph to D. odonohuei.
The ringed morph, in contrast, was found only on the Pacific coast of North America, from Barkley Sound to San Diego (Fig. 1). In the northern portion of its range it was found primarily in subtidal habitats; of the 161 D. sandiegensis found in coastal subtidal habitats 99% were ringed morphs. In addition, all of the 161 D. sandiegensis found south of Fort Bragg, California, were ringed morphs. In most of California, where the spotted morph was absent, the ringed form occupied intertidal, subtidal and bay habitats (Fig. 3). The ringed morphs seemed gradually to decrease in frequency in the intertidal zone north of Trinidad, California (Fig. 3).
Sponge use
Field investigations revealed that both morphs from coastal intertidal habitats fed primarily on the intertidal sponge Haliclona sp. A. Of the 178 individuals found in the intertidal zone in 2009, 51% were found on or within 6 cm of sponge that showed signs of having been fed upon. Of the 91 individuals found on or near sponge, 96% were associated with Haliclona sp. A. Two spotted morphs from the intertidal zone of Trinidad, California, appeared to have fed on a yellow sponge (not identified) and two ringed morphs from Yaquina, Oregon, appeared to have fed on a nearly white sponge (not identified). At intertidal sites, frequencies of the two morphs on or near Haliclona sp. A did not differ (χ² = 2.289, df = 1, P = 0.130).
Faecal analysis indicated that the ringed morph from the subtidal of Monterey Bay fed on Neopetrosia problematica. The primary spicule found in the faeces of seven nudibranchs was a strongyle with mucronate ends, which corresponded to the primary spicule of N. problematica, whereas none of the spicules found in their faeces corresponded to those of Haliclona sp. A.
Reproductive compatibility
Mating studies in the laboratory revealed positive assortative mating among morphs. During two mating trials, individuals that encountered one another for the first time always mated with like morphs (n = 11 matings) (Table 2). In 319 total hours of time-lapse photography nudibranchs continued to choose like morphs in 53 out of 54 repeat mating events, in which 12 out of the 19 nudibranchs mated with multiple partners, and 11 of the 19 mated more than once with the same individuals (Table 2). Mating events lasted from 9 min to 35.5 h (mean 9 h). There were no observable premating behaviours. Each nudibranch circled the aquarium, climbing over multiple other individuals of both morphs, until finding a mating partner. In the no-choice run with four isolated pairs of dissimilar morphs, 22.3 h of time-lapse photography showed only one mating event, lasting for 2.8 h.
Mating events between morphotypes of Diaulula sandiegensis s. l. at field sites from northern California and Oregon where ranges of morphotypes overlap (2009, 2010, 2011) and matings during two mating trials (trial 1 in 2010, trial 2 in 2011) in laboratory between nudibranchs that had met for the first time (initial trials) or had repeatedly mated (extended runs using same individuals as initial trials).
| Field observations | Laboratory trials | |||||
|---|---|---|---|---|---|---|
| Trial 1 (2010) | Extended runs from Trial 1 (87 h) | Trial 2 (2011) | Extended runs from Trial 2 (232 h) | Total matings | ||
| Number of individuals | ||||||
| Spotted | 262 | 6 | 3–6 per run | 4 | 4 per run | |
| Ringed | 46 | 4 | 3–4 per run | 4 | 4–5 per run | |
| Number of unique mating events | ||||||
| Spotted–spotted | 25 (96%) | 2 | 3 | 5 | ||
| Ringed–ringed | 1 (4%) | 3 | 3 | 6 | ||
| Spotted–ringed | 0 (0%) | 0 | 0 | 0 | ||
| Number of repeat mating events | ||||||
| Spotted–spotted | 0 | 26 | 26 | |||
| Ringed–ringed | 2 | 25 | 27 | |||
| Spotted–ringed | 0 | 1 | 1 | |||
| Field observations | Laboratory trials | |||||
|---|---|---|---|---|---|---|
| Trial 1 (2010) | Extended runs from Trial 1 (87 h) | Trial 2 (2011) | Extended runs from Trial 2 (232 h) | Total matings | ||
| Number of individuals | ||||||
| Spotted | 262 | 6 | 3–6 per run | 4 | 4 per run | |
| Ringed | 46 | 4 | 3–4 per run | 4 | 4–5 per run | |
| Number of unique mating events | ||||||
| Spotted–spotted | 25 (96%) | 2 | 3 | 5 | ||
| Ringed–ringed | 1 (4%) | 3 | 3 | 6 | ||
| Spotted–ringed | 0 (0%) | 0 | 0 | 0 | ||
| Number of repeat mating events | ||||||
| Spotted–spotted | 0 | 26 | 26 | |||
| Ringed–ringed | 2 | 25 | 27 | |||
| Spotted–ringed | 0 | 1 | 1 | |||
Mating events between morphotypes of Diaulula sandiegensis s. l. at field sites from northern California and Oregon where ranges of morphotypes overlap (2009, 2010, 2011) and matings during two mating trials (trial 1 in 2010, trial 2 in 2011) in laboratory between nudibranchs that had met for the first time (initial trials) or had repeatedly mated (extended runs using same individuals as initial trials).
| Field observations | Laboratory trials | |||||
|---|---|---|---|---|---|---|
| Trial 1 (2010) | Extended runs from Trial 1 (87 h) | Trial 2 (2011) | Extended runs from Trial 2 (232 h) | Total matings | ||
| Number of individuals | ||||||
| Spotted | 262 | 6 | 3–6 per run | 4 | 4 per run | |
| Ringed | 46 | 4 | 3–4 per run | 4 | 4–5 per run | |
| Number of unique mating events | ||||||
| Spotted–spotted | 25 (96%) | 2 | 3 | 5 | ||
| Ringed–ringed | 1 (4%) | 3 | 3 | 6 | ||
| Spotted–ringed | 0 (0%) | 0 | 0 | 0 | ||
| Number of repeat mating events | ||||||
| Spotted–spotted | 0 | 26 | 26 | |||
| Ringed–ringed | 2 | 25 | 27 | |||
| Spotted–ringed | 0 | 1 | 1 | |||
| Field observations | Laboratory trials | |||||
|---|---|---|---|---|---|---|
| Trial 1 (2010) | Extended runs from Trial 1 (87 h) | Trial 2 (2011) | Extended runs from Trial 2 (232 h) | Total matings | ||
| Number of individuals | ||||||
| Spotted | 262 | 6 | 3–6 per run | 4 | 4 per run | |
| Ringed | 46 | 4 | 3–4 per run | 4 | 4–5 per run | |
| Number of unique mating events | ||||||
| Spotted–spotted | 25 (96%) | 2 | 3 | 5 | ||
| Ringed–ringed | 1 (4%) | 3 | 3 | 6 | ||
| Spotted–ringed | 0 (0%) | 0 | 0 | 0 | ||
| Number of repeat mating events | ||||||
| Spotted–spotted | 0 | 26 | 26 | |||
| Ringed–ringed | 2 | 25 | 27 | |||
| Spotted–ringed | 0 | 1 | 1 | |||
At field sites where the ranges of spotted and ringed morphs overlapped, 26 mating events between unique individuals were photographed. Of these, 25 were spotted with spotted and 1 was ringed with ringed (Table 2). At Yaquina Point, Oregon, there were more of the spotted than of the ringed morph (n = 32 and 2, respectively), but nevertheless the two ringed morphs were mating with one another.
Phylogenetic analyses and haplotype network
Bayesian phylogenetic tree of the Diaulula sandiegensis species complex with D. hummelincki as the outgroup, based on concatenated sequences of mitochondrial COI and 16S genes. Posterior probabilities (PP > 0.5) are shown above the branches and maximum-likelihood bootstrap (BS > 50) values shown below. Ringed morph corresponds to D. sandiegensis s. s., spotted morph to D. odonohuei. Membership of the haplogroups of D. odonohuei recovered in the haplotype network analysis (Fig. 5) is indicated. Specimens coded with symbols based on geographic location; black triangles: eastern Pacific, black circles: western Pacific.
Haplotype network for COI sequence data of Diaulula sandiegensis species complex showing two distinct networks. A. Haplotype network for D. sandiegensis s. s. B. Haplotype for D. odonohuei with three distinct groups of haplotypes (A1–A3). Each circle represents a unique haplotype and its size indicates the number of specimens sharing that haplotype. Each line between circles indicates a single substitution; small white circles represent hypothetical haplotypes. Ancestral haplotypes marked with an asterisk.
ABGD analysis
Using K2P, TN and uncorrected distance matrices, the COI and 16S sequences showed a barcode gap between a priori genetic distance thresholds of 0.01 and 0.02. Using a value of P between this range (0.0129), two species groups were identified. For COI, assignment of individuals to the two groups matched the two colour-morph clades from the Bayesian and ML analyses, but not for 16S. However, using a slightly higher value (0.0215), the same two groups were identified for 16S as well. The geographic distribution of the two species is shown in Figure 1.
Molecular-clock analysis
Divergence times from the molecular-clock analyses are summarized in Table 3. The Langley–Fitch method, which estimates one substitution rate across the entire tree, yielded a rate of 0.25% per myr for COI, which is similar to other fossil-calibrated phylogenies of tropical gastropods (Frey & Vermeij, 2008; Reid, Dyal & Williams, 2010). Two nodes were dated, the first (1.5–1.76 Ma) between the two species here recognized and the second (0.31–0.35 Ma) between one of the western Pacific clades of the transpacific species (haplogroup A1) and the clade containing both eastern and western Pacific specimens (haplogroups A2 + A3) (see Supplementary Material Fig. S5).
Results from r8s molecular-clock analyses, including mean and standard deviation (SD) for divergences (millions of years ago, Ma) using different methods and algorithms. For explanation of node codes see Fig. 6.
| Node | Langley–Fitch/Local molecular clock | Langley–Fitch/Powell | Langley–Fitch/quasi-Newton | Langley–Fitch/truncated Newton | Non parametric rate smoothing/Powell | Penalized likelihood/Powell | Penalized likelihood/truncated Newton |
|---|---|---|---|---|---|---|---|
| A | 1.5 | 1.5 | 1.5 | 1.5 | 1.76 | 1.76 | 1.76 |
| B | 0.31 | 0.31 | 0.31 | 0.31 | 0.34 | 0.35 | 0.35 |
| SD | ±1.03 × 10−8 | ±2.05 × 10−8 | ±1.37 × 10−8 | ±1.37 × 10−8 | ±0.278 | ±0.00922 | ±0.00922 |
| Node | Langley–Fitch/Local molecular clock | Langley–Fitch/Powell | Langley–Fitch/quasi-Newton | Langley–Fitch/truncated Newton | Non parametric rate smoothing/Powell | Penalized likelihood/Powell | Penalized likelihood/truncated Newton |
|---|---|---|---|---|---|---|---|
| A | 1.5 | 1.5 | 1.5 | 1.5 | 1.76 | 1.76 | 1.76 |
| B | 0.31 | 0.31 | 0.31 | 0.31 | 0.34 | 0.35 | 0.35 |
| SD | ±1.03 × 10−8 | ±2.05 × 10−8 | ±1.37 × 10−8 | ±1.37 × 10−8 | ±0.278 | ±0.00922 | ±0.00922 |
Results from r8s molecular-clock analyses, including mean and standard deviation (SD) for divergences (millions of years ago, Ma) using different methods and algorithms. For explanation of node codes see Fig. 6.
| Node | Langley–Fitch/Local molecular clock | Langley–Fitch/Powell | Langley–Fitch/quasi-Newton | Langley–Fitch/truncated Newton | Non parametric rate smoothing/Powell | Penalized likelihood/Powell | Penalized likelihood/truncated Newton |
|---|---|---|---|---|---|---|---|
| A | 1.5 | 1.5 | 1.5 | 1.5 | 1.76 | 1.76 | 1.76 |
| B | 0.31 | 0.31 | 0.31 | 0.31 | 0.34 | 0.35 | 0.35 |
| SD | ±1.03 × 10−8 | ±2.05 × 10−8 | ±1.37 × 10−8 | ±1.37 × 10−8 | ±0.278 | ±0.00922 | ±0.00922 |
| Node | Langley–Fitch/Local molecular clock | Langley–Fitch/Powell | Langley–Fitch/quasi-Newton | Langley–Fitch/truncated Newton | Non parametric rate smoothing/Powell | Penalized likelihood/Powell | Penalized likelihood/truncated Newton |
|---|---|---|---|---|---|---|---|
| A | 1.5 | 1.5 | 1.5 | 1.5 | 1.76 | 1.76 | 1.76 |
| B | 0.31 | 0.31 | 0.31 | 0.31 | 0.34 | 0.35 | 0.35 |
| SD | ±1.03 × 10−8 | ±2.05 × 10−8 | ±1.37 × 10−8 | ±1.37 × 10−8 | ±0.278 | ±0.00922 | ±0.00922 |
DISCUSSION
Taxonomy
Based on the molecular data obtained in this study, the name Diaulula sandiegensis, as currently used in the modern literature (e.g. Behrens & Valdés, 2001), represents a complex of two ‘pseudocryptic’ species (distinguishable morphologically but only following molecular distinction; Hoover et al., 2015). This conclusion was initially supported by the phylogenetic and ABGD analyses of mitochondrial sequence data, but after a reevaluation of morphological traits the external colour differences between the two species became obvious. The ‘spotted morph’ species exhibits a large number of irregular chocolate-brown spots that can be ring shaped or solid and extend onto the mantle margin (Fig. 2A, B, E–G). In contrast, the ‘ringed morph’ species has a small number of large chocolate-brown ring-shaped spots located only on the dorsum and spotting does not extend onto the mantle margin (Fig. 2C, D, H). Behavioural observations suggest that the species are reproductively isolated, as morph-specific mate choice was the dominant pattern in mating trials.
Doris sandiegensis was originally described by Cooper (1863), based on specimens collected in San Diego Bay, California, as a pale brownish yellow species with 12–20 annular or entirely brown spots. Although the description is brief and did not include illustrations of the live animals, this colour pattern is clearly distinct from those of other northeastern Pacific dorids. Therefore, this name has been consistently and correctly applied to the species (McDonald, 1983), at least in the southern portion of the range. Later, the species was used as the type of the new genus Diaulula by Bergh (1878). The validity of Diaulula (a member of the caryophyllidia-bearing dorid clade) was later confirmed by Valdés & Gosliner (2001), based on a morphological phylogenetic analysis. Because only the ringed morph occurs in southern California, including the type locality of D. sandiegensis, this name is retained herein for the southern species.
A review of the literature reveals that there is a possible available name for the northern, transpacific species, D. odonoghuei. O'Donoghue (1922: 149) described Doris echinata from Vancouver Island as an opaque white species with 12–40 small spots of a “warm brown colour” scattered irregularly over the dorsal surface. O‘Donoghue (1922) also described the caryophyllidia in this species as blunt conical papillae with bundles of spicules passing up the sides of each papilla, as well as the radula and other anatomical details, all of which are similar to the characteristics of D. sandiegensis. Subsequently, O'Donoghue (1926: 206) pointed out that “Doris echinata is really Doridigitata maculata” (a new species name) without further explanation. Steinberg (1963) indicated that both Doris echinata O'Donoghue, 1922 and Doridigitata maculata O'Donoghue, 1926 were preoccupied, by Doris echinata Lovén, 1846 and Doris maculata Garstang, 1896, respectively; and therefore, she proposed the replacement name Doris odonoghuei for this species. Behrens & Valdés (2001) considered Doris odonoghuei Steinberg, 1963 as a synonym of D. sandiegensis. Attempts to locate the type material of Doris odonoghuei failed; the specimens are not housed at the Royal Ontario Museum, the Canadian Museum of Nature, Ottawa or the Royal British Columbia Museum (formerly Provincial Museum of Victoria). Specimens examined from Barkley Sound (near the type locality of Doris odonoghuei, Vancouver Island) are of both morphs, but in this region the spotted morph is the only one present in intertidal environments, which is where O'Donoghue would have collected specimens. O'Donoghue's (1922) description of 12–40 small brown spots scattered irregularly over the dorsal surface for D. echinata is inconclusive, but more consistent with the characteristics of the spotted morph than with those of the ringed morph. The identity of Doris odonoghuei is thus ambiguous. Therefore, and to avoid the introduction of a new name for the spotted form, we propose to fix the use of the name Diaulula odonoghuei for the northern, transpacific species by the designation of a neotype. The neotype here designated (LACM 3349) was collected in Gig Harbor, Washington, and its identity as a member of the northern, transpacific species was verified with sequence data (Supplementary Material Table S1: isolate TL027).
Biogeography
Diaulula odonoghuei has a large geographic range across the North Pacific, from Korea (Koh, 2006) and Onagawa, Japan (Baba, 1957) to the Pacific coast of North America southwards to Bodega Bay in northern California (Kelly, 2013). In contrast, D. sandiegensis has a more restricted range along the Pacific coast of North America, from the outer coast of Baja California (Behrens & Hermosillo, 2005) at least to Barkley Sound, British Columbia (Fig. 3). The ranges of the two species overlap between Fort Bragg and British Columbia, although D. sandiegensis is difficult to find from northern California to Washington, only to appear again subtidally in Tatoosh Island, Washington, Barkley Sound and other areas in British Columbia (Penney, 2013a; Fig. 3).
Diaulula odonoghuei has an amphi-Pacific range, but because of the absence of shared haplotypes between western (Russia and Japan) and eastern Pacific (North America) localities, it seems unlikely that there is current gene flow between the sampled populations, i.e. by larval dispersal or rafting of adults on floating debris. Similar disjunction has been shown for other benthic marine organisms (Cox, Zaslavskaya & Marko, 2014). Thus, the amphi-Pacific range of D. odonoghuei seems to be the result of past events, either dispersal across the North Pacific or vicariance of a former amphi-Pacific range. The haplotype network analysis (Fig. 5) also recovered a group of haplotypes from Russia and Japan (haplogroup A3) that are substantially divergent from others in the western populations of D. odonoghuei (haplogroup A1), by more than 11 substitutions, grouping instead with the eastern Pacific haplotypes (haplogroup A2). This suggests a complex biogeographic history including multiple dispersal or vicariant events. Further research, using larger sample sizes and more variable markers, may reveal additional cryptic diversity within D. odonoghuei.
The molecular-clock analyses estimate that D. sandiegensis and D. odonoghuei diverged between 1.5 and 1.76 Ma. This timing coincides with a cooling period during the Calabrian age of the Pleistocene (Cita et al., 2008; Gibbard, Head & Walker, 2009). The cold period at around 1.5 Ma (Hays, Imbrie & Shackleton, 1976; Wunsch, 2004; Lisiecki & Raymo, 2005; Laskar et al., 2011) led to the expansion of ice sheets in the North Pacific and could have resulted in an allopatric speciation event between D. sandiegensis and D. odonoghuei. This hypothesis requires that D. sandiegensis evolved in the eastern Pacific and D. odonoghuei in the western Pacific, and that the latter subsequently dispersed eastward. The molecular-clock analyses also indicate that the clade including North American and Sea of Japan specimens (haplogroups A2+A3) of D. odonoghuei separated from that found solely in the Sea of Japan (haplogroup A1) around 0.31–0.35 Ma, during the Ionian age of the Pleistocene. This timing coincides with another cooling period, immediately following a warm period with an average temperature +2 °C above the 1960–1990 average (Jouzel et al., 2007; Hansen et al., 2013). Thus, we hypothesize that the current population structure of D. odonoghuei is the result of dispersal from west to east during one or more interglacial periods followed by vicariance during glaciations. Although the split between haplogroups A2 and A3 could not be dated as it is much too recent, it suggests yet another east-west divergence, providing additional support for the proposed scenario. Additional research including other groups of organisms with similar ecologies is necessary to test and refine the proposed scenario; unfortunately, the absence of a fossil record for nudibranchs will hamper our ability to provide definite support for any biogeographic hypothesis.
Mode of speciation
The fact that the main divergence events between and within D. sandiegensis and D. odonoghuei appear to coincide with past glacial events, potentially resulting in vicariance between east and west, supports the hypothesis that the recent evolution of this group has been shaped by climate change and that the two species diverged allopatrically. Mating trials in the laboratory revealed strong evidence for positive assortative mating by species, which suggests the existence of pre-zygotic isolation. Although no differences in premating behaviours were seen between the two species, individuals that encountered one another for the first time in laboratory trials always mated with the same species (11 matings in trials 1 and 2). When all 92 matings observed in the field and laboratory are combined, 98% occurred within species, even after crawling all over the other species in the laboratory (Table 3). This suggests that a premating barrier, such as chemical recognition, may affect mate choice and more importantly that the speciation process is complete. Although the two species are partially sympatric, D. odonoghuei is restricted to intertidal and bay habitats from Fort Bragg, California, northwards, whereas D. sandiegensis is restricted to subtidal habitats where its range overlaps with that of D. odonoghuei, suggesting character displacement after secondary contact.
Further research on other groups of benthic organisms is necessary to explore whether the Northern Hemisphere glaciations have acted as engines of evolution, resulting in the production of diversity. This has important implications for understanding the future of the North Pacific biota as anthropogenic climate change may lead to further transpacific dispersal.
SUPPLEMENTARY MATERIAL
Supplementary material is available at Journal of Molluscan Studies online.
ACKNOWLEDGEMENTS
Financial support was provided by a California State Polytechnic University Provost Teacher–Scholar Award to AV, an Ernest Prete, Jr. Student Environmental Research Fellowship and a Conchologists of America student award both to TL. TL and HK thank H. Fortunato and L. Eduardo (Hokkaido University) for their support in arranging fieldwork in Japan. JK thanks Mike Kelly and Jen Kelly for their long hours spent searching for nudibranchs in tidepools and Niusha Taeidi-Mackie (San José State University) for assistance with laboratory work. JK and SFC thank M. Fletcher (Humboldt State University), E. Thomson and A. Carter (NSF Research Experience for Undergraduates) for their hard work searching for D. sandiegensis, G. Eberle and D. Hoskins (HSU Telonicher Marine Laboratory) for their support with the setup of saltwater aquaria and A. Baker (Humboldt State University‘s common core genetic facility) for his assistance with genetic analysis. Discussions with E. C. Metz, J. O. Reiss, E. S. Jules, R. Alvarez, S. Monk, K. Korcheck, R. Koeppel and J. Koeppel from Humboldt State University improved this manuscript. Thanks also to B. Grasse from the Monterey Bay Aquarium for his help with collection of nudibranchs from Monterey Bay and both D. Behrens and B. Penney for the initial idea for this study. M. Frey, Royal British Columbia Museum, D. Currie and M. Zubowski, Royal Ontario Museum, and J.-M. Gagnon, Canadian Museum of Nature, assisted with attempts to locate O'Donoghue‘s type material, and L. Groves (LACM) helped with curation of specimens and access to the collection. Two anonymous reviewers as well as Editor D. Reid and Associate Editor P. Marko provided valuable comments that greatly improved the manuscript. This project was supported by grants from the California State University Program for Education and Research in Biotechnology (CSUPERB) to AV, the National Science Foundation (NSF: DBI-0755466) to SFC and JM, and the Far East Research Program (FERP: 15-I-6-014o) to AC. The SEM work was conducted with the assistance of G.-A. Kung at the LACM SEM laboratory supported by the NSF: DBI-0424911 to AV et al.
REFERENCES
Author notes
These two authors contributed equally.
