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).