Materials and methods
Snails were collected in January-March 2015 from coral colonies at ten localities along the leeward coast of Curaçao, southern Caribbean (Fig. 1). Colonies of scleractinians and alcyonaceans were haphazardly selected and searched for snails during SCUBA dives. As the sampling effort was not equal among localities and localities were < 50 km apart from each other, locality was not used as a factor in statistical analyses. Aforementioned surveys of invertebrates associated with corals appear to be effective, as they have previously resulted in new host records for Curaçao and the Caribbean in for example ovulid snails on octocorals (Reijnen et al., 2010), gall crabs in scleractinians (Van der Meij, 2014) and serpulids in reef corals (Hoeksema et al., 2015; Hoeksema and Ten Hove, 2017).
All snails associated with a colony were collected and stored in one plastic sampling bag per host colony. Host corals were photographed and their identity, depth and locality were recorded. If present, feeding scars were photographed as well. After sampling, snails were put in 96% ethanol (in a few cases 70% ethanol) until further processing. All specimens have been deposited in the collections of Naturalis Biodiversity Center, Leiden, the Netherlands (catalogued under numbers coded RMNH). In the case of snails associated with alcyonaceans, a small sample of the distal part of host colonies was also collected for species identification. Alcyonacean sclerites were isolated by dissolving the coenchymal tissue in sodium hypochlorite (4% household bleach solution). Alcyonaceans were identified to the genus level using photographs and light microscopy slides of the sclerites, using the keys by Bayer (1961).
Shell length of the snails was measured with a digital calliper to the nearest 0.01 mm. Of 60 snails (7% of total), only the shell length was measured with a Vernier calliper to the nearest 0.05 mm. Shell length was defined as the length from the tip of the apex to the tip of the aperture (end of the anterior canal) (Fig. 2). To determine any measurement error that could arise from inconsistencies in the orientation of the shell between the calliper blades, the shell length of 65 snails was measured in triplicate. The measurement error was defined as the average distance to the mean of the three replicate measurements of each shell. The error for shells measured with the Vernier calliper was not calculated. Shell lengths and widths were log-transformed in all statistical analysis to achieve normal distributions and homogeneity of variance. Differences in shell length were tested using ANOVA models and in some cases using Kruskall-Wallis rank sum tests.
Fig. 2. Shell length as used to measure all shells. Twelve black points represent the 12 landmarks used in the geometric morphometric analysis.
Landmark-based geometric morphometrics were used to assess the shape of each shell. After measuring shell length, the ventral side (aperture facing upwards) was photographed with a Nikon D7000 DSLR camera equipped with a Sigma 105 mm macro lens. The locations of 12 landmarks were recorded on each photo (Fig. 2) and chosen to capture the observed variation in shell shape during shell measurements. Most of these landmarks have been used before in the morphometrics of gastropods (Zelditch et al., 2004; Hollander et al., 2006; Mariani et al., 2012). Shells covered by encrusting algae were excluded from this analysis because their landmarks were hidden. To align landmarks and remove the effect of size, a generalized Procrustes superimposition was applied to the data (Gower, 1975; Rohlf and Slice, 1990). Replication errors landmark data were calculated as well. See Online Supplementary Material 1 for the methods followed.
To statistically test for differences in shell shape and potential relationships between shell shape and shell size, snail host species or depth (as well as the interactions between host species and both shell size and depth), distance-based Procrustes ANOVA models were used that are equivalent to other distance-based ANOVA methods, like PerMANOVA (Goodall 1991; Anderson 2001; Adams and Otárola-Castillo 2013). For all Procrustes ANOVA models, significance of the different factors was tested against 10,000 permutations. Host-associated differences in shell shape were tested through pairwise comparisons of the effect of host species in a full model (with all tested factors included) against a reduced model (with all factors except host species included). To account for multiple tests, p-values were corrected using a Bonferroni correction.
To define allometric patterns, a common allometric component (CAC) was calculated from the landmark data to express allometric patterns as one variable (Mitteroecker et al., 2004). Host-specific regressions between CAC and shell length were made (excluding hosts having less than five specimens with morphometric data). The vectors of shell length of snails associated with different host species were compared to reveal differences in allometric patterns in the amount of change in shell shape per unit of growth (corresponding to the distance among vectors of shell length) and the direction of shell shape change (corresponding to the correlation among vectors of shell length). Pairwise comparisons of both the distance and correlation among the vectors of shell length were made between a full model and a reduced model without the interaction between the factors host species and shell length (see Online Supplementary Material 2). As before, the p-values were corrected with a Bonferroni correction.
To visualize variation in shell shape, a principal component analysis (PCA) was performed using the landmark data. To separate real variation in shell shape from noise resulting from the error described above and calculate repeatability of axes, landmark data of the three replicated photos was included in the PCA. Intraclass correlation coefficients (ICC, model 2,1) were calculated between the PCA-scores of triplicates on all PCA-axes. An axis was considered repeatable when the ICC was higher than 0.80 (Burridge et al., 2015). To visualize differences in allometric patterns among snails associated with different host species, linear regressions between PCA scores and shell length were used to predict PCA scores (and therefore shell shape) of shells of specific lengths associated with specific host species
After the morphometric analysis, shells were crushed to remove the snail from its shell. Using a dissecting microscope, the sex of each snail was determined by presence or absence of a penis just above the left eyestalk. Differences in sex ratios among snails associated with different host species were assessed using Fisher’s exact tests. For tests on larger tables (to test for differences among host species), p-values were computed using Monte Carlo simulations, with 1 million replications. Pairwise differences among host species were assessed using pairwise Fisher’s exact tests; p-values were corrected using a Bonferroni correction. Linear regressions were made between sex ratio and mean shell length, and between mean male and female shell length. Individual points (corresponding to a single host species) were weighed according to the number of specimens. Only host species with more than five specimens were included in the analyses comparing snails among host species.
A small piece of tissue was removed from the foot of a selection of snails for genetic analysis. Two mitochondrial markers were amplified and sequenced: a fragment the 12S rRNA (12S) gene and a fragment of the cytochrome c oxidase subunit 1 gene (COI). Both markers have been used extensively and proven informative in closely related gastropods (Oliverio and Mariottini, 2001a, 2001b; Barco et al., 2010; Gittenberger and Gittenberger, 2011). DNA was extracted on a KingFisher Flex magnetic particle processor (Thermo Scientific), using the Nucleospin Tissue kit (Macherey-Nagel, Düren, Germany). PCR reaction mixtures consisted for both markers of 0.25 µL QIAGEN Taq DNA polymerase (5 units µL-1), 0.5 µL dNTPs (2.5 mM) and 1.0 µL of both the forward and reverse primers, as well as 0.5 µL 100 mM Promega BSA, 0.5 µL 25 mM MgCl2 , 2.5 µL 10x PCR buffer (QIAGEN) and 15.8 µL milli-Q water. In the PCR, an annealing temperature of 50°C was used for both 12S and COI. PCR products were sequenced using Sanger sequencing by BaseClear (Leiden, the Netherlands). Primers were used as in Barco et al. (2010) (Table 1). Five previously published sequences were used in the phylogenetic analysis (Table 2).
Table 1. Sequences of forward (F) and reverse (R) primers for the amplification of two mitochondrial markers, 12S rRNA (12S) and cytochrome c oxidase subunit I (COI). All primers are as in Barco et al. (2010).
Table 2. Specimens used in phylogenetic analysis with their voucher numbers (RMNH) and GenBank accession numbers.
Forward and reverse sequences were assembled automatically and edited by hand in Sequencher 5.4 (Gene Codes Corporation, Ann Arbor, MI, USA). Sequences were aligned using the MAFFT algorithm on the GUIDANCE2 server (Katoh et al., 2005; Sela et al., 2015). The appropriate substitution model for either marker was determined based on the Akaike information criterion (AIC) in both the jModelTest 2.1.7 and MrModeltest 2 software packages (Nylander, 2004; Darriba et al., 2012). Both software packages agreed on the best substitution model to be used in MrBayes. The GTR + Γ (Tavaré 1986) model was used for 12S, the HKY85 + Γ + I (Hasegawa et al., 1985) model was used for COI. A phylogenetic tree was constructed based on both markers separately as well as a concatenated dataset using Bayesian inference with the (parallel) Metropolis coupled Markov chain Monte Carlo ((MC)3) method in MrBayes 3.2 (Altekar et al., 2004; Ronquist et al., 2012). In MrBayes, the (MC)3 analysis was run in duplicate, for a length of 25,000,000 generations for the concatenated dataset and a length of 15,000,000 generation for the trees based on single markers. Trees were sampled every 100 generations. Burn-in was determined by looking at the deviation of split frequencies, trees sampled before this deviation dropped below 0.01 were discarded. For the concatenated analysis, the burn-in was determined to be 1.49 million generations, almost 6% of the total length of the analysis. For the trees based on 12S and COI separately, burn-in was determined to be 555,000 and 1,155,000 generations respectively. Sequences of the muricid species Drupella rugosa (Born, 1778), previously published by Claremont et al. (2011), were used as an outgroup.
Finally, an automatic barcode gap discovery analysis (ABGD) based on Kimura two-parameter (K2P) model on the marker COI was used to assess species delineation on the phylogenetic tree and to identify Molecular Operational Taxonomic Units (MOTUs) (see Kimura, 1980; Blaxter, 2004; Puillandre et al., 2012; Barco et al., 2013). Sequences from coralliophiline snails previously published by Harasewych et al. (1997), Puillandre et al. (2009), Barco et al. (2010), Claremont et al. (2011) and Gittenberger and Gittenberger (2011) were included in this analysis (Online Supplementary Material 3), while the outgroup was excluded.
To further assess host-associated genetic divergence, a haplotype network was built for both C. galea and C. caribaea, using both markers separately. Networks were calculated using an infinite site model based on uncorrected distances between haplotypes. To statistically assess intraspecific genetic divergence, an AMOVA (p-values calculated based on 100,000 permutations) with host species or host order for both C. galea and C. caribaea was used on both markers. Statistics on genetic data and the calculation of the haplotype networks were performed in R using the packages APE 3.4 and pegas 0.8-2 (Paradis et al., 2004; Paradis, 2010; R Core Team, 2015).