Population structure of the fireworm Hermodice carunculata in the wider Caribbean, Atlantic and Mediterranean Sea

Abstract

The bearded fireworm Hermodice carunculata is a conspicuous omnivorous benthic invertebrate distributed in the wider Atlantic Ocean, Mediterranean Sea and the Red Sea. Questions were raised in the past regarding the taxonomic status of amphiatlantic fireworm specimens but subsequent analysis of morphological data with new molecular data reconfirmed the presence of only one species, H. carunculata with amphiatlantic distribution. Hermodice carunculata has been the subject of several taxonomic and genetic studies but there are still questions regarding the genetic diversity and population structure of the species throughout its range. We contributed to the genetic studies of H. carunculata with ~200 new samples from over 20 locations in the Caribbean, Eastern Atlantic Ocean and the Mediterranean Sea. We sequenced the mitochondrial markers Cytochrome c oxidase Subunit I (COI) and Cytochrome b (Cytb) from each polychaete to examine the patterns of genetic diversity of H. carunculata. Our data revealed a significant population structure between the Caribbean H. carunculata including those from Brazil compared with those from the eastern Atlantic (Tenerife, Canary Islands) and eastern Mediterranean (Greece and Malta). We inferred rates of gene flow between eastern and western Atlantic populations by estimating the bidirectional effective migration rate (N_em). We corroborate previous studies indicating the existence of one species with genetically distinct populations. However, the biological significance of the observed population divergence should be evaluated with cross-breeding experiments.

Introduction

Amphinomid polychaetes are among the most recognized epibenthic polychaetes because of the burning sensation caused by their chaetae after penetrating human skin (Yáñez-Rivera & Salazar-Vallejo, 2011); thus, their common name ‘fireworms’. They are found in tropical and subtropical waters with a relatively wide depth distribution (0–300 m), but they are encountered mostly in the intertidal zone (Pardo & Amaral, 2006; Yáñez-Rivera & Salazar -Vallejo, 2011; Wolf & Nugues, 2013). Many species are cryptic as they are abundant in coral reefs or rocky areas, where crevices provide them with suitable habitats (Pardo & Amaral, 2006). They are also found in mudflats and turtle grass beds (Sharp et al., 2005).

In the wider Caribbean, the most well-known fireworm species is Hermodice carunculata (Pallas, 1766), also called the bearded fireworm. Originally described as Aphrodita carunculata based on the shape and development of the olfactory organ caruncle, it was later moved to the genus Hermodice by Kinberg in 1857 (Yáñez-Rivera & Salazar-Vallejo, 2011). Hermodice carunculata displays a high degree of dietary plasticity from detritus to live or injured/dead metazoans such as sponges, zoanthids, scleractinian corals, sea anemones, gorgonians, milleporid hydrocorals, sea urchins, starfish, sea cucumbers and fish (Witman, 1988; Sharp et al., 2005; Pardo & Amaral, 2006; Pérez & Gomes, 2012; Barroso et al., 2016, 2017; Schulze et al., 2017).

The fireworm's efficient way of scavenging was probably not of a concern when reefs were still pristine but nowadays, fireworms may lead to further declines in the depauperate coral reefs of the Caribbean (Wolf et al., 2014; Schulze et al., 2017). Past studies have shown how H. carunculata can cause permanent damage, slow recovery and vertical growth of corals, as well as promote algal colonization, thus aiding in the baseline shift of tropical reefs (Hughes, 1994). Witman (1988) estimated that H. carunculata open 12.9 cm² of new space per 1.0 m² on the hydrocoral Millepora complanata per day. Only 22–25% of the lesion space was recovered by tissue growth after one year and the vertical growth of Millepora branch tips predated by H. carunculata was significantly lower than non-predated branch tips (Witman, 1988). Furthermore, H. carunculata reinforces recruitment bottlenecks and reduces the recovery from disturbance, since they particularly feed on new reef settlers causing early post-settlement mortality in corals (Wolf & Nugues, 2013). For instance, the predation on both adults and newly settled corals threatens the existence and the potential recovery of species such as the reef-building Orbicella spp. (Wolf & Nugues, 2013). The bearded fireworms have been suggested as disease vectors and pathogen transmitters, increasing vulnerability to coral disease and/or bleaching as they act as a winter reservoir and summer vector for the pathogen Vibrio shiloi in the Mediterranean (Sussman et al., 2003). In some cases, the feeding scars of H. carunculata on corals have been interpreted as coral disease and the ecological impact of the fireworms may have been underestimated (Miller et al., 2014). There are few predators known to prey on H. carunculata; e.g. some Caribbean cone snails (Vink, 1974; Vink & von Cosel, 1985), white grunts and sand tilefish (Ladd & Shantz, 2016). There have also been observations of bluehead (Thalassoma bifasciatum) and yellowhead wrasses (Halichoeres harnoti) preying on H. carunculata (Wolf, 2012) but most likely, there is little or no effect of predation on the abundance of fireworms.

Hermodice carunculata is considered the sole representative of the genus but the global phylogeography of the species has not been thoroughly examined. Initially, Baird (1870) recognized six species in Hermodice and described H. nigrolineata from the Mediterranean Sea and the Canary Islands, which was later synonymized with H. carunculata, as was H. savignyi (Fauvel, 1923; Ebbs, 1966). Hermodice nigrolineata was re-established when H. carunculata was redescribed from the topotype specimen by Yáñez-Rivera & Salazar-Vallejo (2011). According to Yáñez-Rivera & Salazar-Vallejo (2011), H. carunculata is restricted to the Grand Caribbean region and H. nigrolineata is distributed in the Mediterranean Sea and adjacent eastern Atlantic Ocean areas. A re-examination of a key character (i.e. branchial filaments) used to differentiate H. nigrolineata from H. carunculata questioned the validity of that character to diagnose the eastern and western Atlantic fireworms as different taxa; also, morphology was discordant with the observed genetic population structure (Ahrens et al., 2013). Two more recent studies have shown that the structure of the branchial filaments, the main distinguishing characteristic between H. carunculata and H. nigrolineata according to Yáñez-Rivera & Salazar-Viejo (2011), is highly plastic and therefore not useful for species delineation (Grimes et al., 2020; Lucey et al., 2020). The hypothesis that the eastern Atlantic specimens of Hermodice belong to a different species has been rejected and once again, H. nigrolineata is not recognized as a distinct species (Ahrens et al., 2013).

Furthermore, a phylogenetic analysis of additional COI sequences of H. carunculata from the eastern Mediterranean revealed two main clades, one exclusively found in the Mediterranean and a second comprised surprisingly by both Mediterranean and eastern Atlantic sequences (Chatzigeorgiou et al., 2014). Based on previous studies and the samples collected during our former expeditions, we aimed at updating the existing knowledge on the species phylogeography. Hermodice carunculata is an ideal species to use for genetic analysis to infer patterns of population structure and gene flow because of its ecological importance, amphiatlantic distribution and its biphasic lifestyle with planktonic larvae (Toso et al., 2020) and benthic, relatively sedentary adults.

Methods

Hermodice carunculata specimens were collected from seven locations in the Caribbean (Curaçao, Dominican Republic, Guadeloupe, Panama, Puerto Rico and St. Croix) eastern Atlantic Ocean (Tenerife of Canary Islands) and eastern Mediterranean Sea (Ikaria Island of Greece) (Table 1, Figure 1). Hereafter, the terms ‘location’ and ‘population’ will be used interchangeably.

Fig. 1. Map with locations denoted with asterisks indicating the geographic spread of samples used in this study. Detailed information on the origin of samples is listed in Table 1.

Table 1. Number of DNA sequences of COI by geographic location used in this study

Asterisks denote locations where we included previously published sequences from Ahrens et al. (2013) and Chatzigeorgiou et al. (2014).

All specimens were hand-collected by scuba or snorkelling in depths ranging from 3 m to 20 m. Even though H. carunculata has been reported mostly active from dusk to dawn (Marsden, 1960), we observed a wide variation in the time of peak movement outside their cryptic habitats during sample collection. Usually, there were very few specimens throughout the day; sometimes they were not seen at night in the same location where they were abundantly found during the day (pers. obs., all in Puerto Rico). There were also significant habitat differences in peak activity among the sampled locations in the Caribbean, Tenerife and Greece. In Tenerife, the fireworms were collected among large barren rocky outcrops in the middle of the day and in Ikaria, Greece, they came out in large numbers after 5:00 pm during the summer months. The Ikarian habitat was 1–2 m deep and was characterized by large boulder rocks covered by macroalgae.

Immediately after collection, all samples were stored in 95% ethanol. DNA was extracted from 1–2 segments of the polychaetes using the Qiagen® Dneasy® Blood and Tissue kit (Qiagen Corp.) according to the manufacturer's instructions. We also used Chelex® 100 Resin (Bio-Rad, Inc.) as an extraction method for 1/4 of the samples. The quantity and quality of DNA were obtained with a NanoDrop 2000™ Spectrophotometer.

We amplified by PCR (Polymerase Chain Reaction) fragments of ~600 bp of the mitochondrial Cytochrome Oxidase subunit I (COI) gene with the universal primers LCOI-1490 (5′ GGT CAA CAA ATC ATA AAG ATA TTG G 3′) and HCOI-2198 (5′ TAA ACT TCA GGG TGA CCA AAA AAT CA 3′) (Folmer et al., 1994). We also amplified a portion of the mitochondrial Cytochrome b gene (Cytb) using the primers Cytb-424F (5′ GGW TAY GTW YTW CCW TGR GGW CAR AT 3′) and Cytb-876R (5′ GCR TAW GCR AAW ARR AAR TAY CAY TCW GG 3′) from Merritt et al. (1998). Each PCR reaction mixture (total volume 25 μl) contained 1 μl DNA (100 ng) as template, 0.3 μl of each primer forward and reverse (100 μM), 12.5 μl of BioMix (Bioline Company) and 10.9 μl water (ddH₂O). Other reaction mixes used were: BioMix Red, 5× MyTaq™ Mix using standard protocols from the manufacturer. PCR reactions were conducted on a Bio-Rad MyCycler ™ Thermal Cycler. The PCR protocol for COI was 95°C for 3 min, followed by 35 cycles of 15 s at 95°C, 30 s at 43°C, 30 s 72°C, followed by 1 cycle of 5 min at 72°C, and finally keeping the PCR products at 15°C. The PCR protocol for Cytb was 95°C for 3 min, followed by 35 cycles of 15 s at 95°C, 30 s at 45°C, 1 min 72°C, followed by 1 cycle of 5 min at 72°C. Products of PCR were then inspected by gel electrophoresis. PCR products were loaded on a 1% agarose gel stained with ethidium bromide and visualized by UV fluorescence. All successful PCR products were processed for sequencing using the Big Dye 3.1 Terminator Cycle Sequencing Kit and the ethanol-precipitated products were loaded for sequencing into an ABI 3130xl 16-capillary Genetic Analyzer at the University of Washington genomic facility (htseq.org). Almost all sequences were sequenced with the forward primer after preliminary sequencing with both primers yielding high-quality DNA traces. In the analysis, we included additional publicly available COI sequences from Ahrens et al. (2013; KC017475-KC017594) and from Chatzigeorgiou et al. (2014; KF878397-KF878476) to increase our geographic coverage. For example, in the Mediterranean COI data set, we have included our own sequences from Greece, as well as those sequences from Crete (Greece) from both Ahrens et al. (2013) and Chatzigeorgiou et al. (2014), and those from Malta (Ahrens et al., 2013).

The resulting chromatograms were viewed and checked for mutations in Codon Code Aligner v. 5.1.5 (Codon Code Corp.). End-trimmed sequences were aligned within Codon Code Aligner using Clustal Ω alignment program. Verified sequences were imported in DNAsp v. 5.10.01 and were grouped by location as listed in Table 1. An Arlequin file .arp (project file) and .hap (file with all haplotypes) were created through DNAsp for downstream analysis. Molecular diversity indices (π, θ), neutrality tests (Tajima's D and Fu's Fs) and population structure (AMOVA test) were estimated in Arlequin 3.5.2.2 (Excoffier & Lischer, 2010). Because of the uneven number of samples per location we estimated the effective number of haplotypes by using the formula n _e = 1/Σp_i ², where p_i is the frequency of the ith haplotype and n _e is the effective number of haplotypes (Crawford et al., 2004). An AMOVA test was conducted between eastern and western Atlantic specimens to test for differentiation across the Atlantic Ocean with our new data and locations. For the AMOVA analysis, the eastern Atlantic samples comprised Tenerife separately from the pooled eastern Mediterranean sites since they are separated by two well-known biogeographic boundaries: the Almeria-Oran front (e.g. Rίos et al., 2002) and the Siculo-Tunisian strait (e.g. Zardoya et al., 2004) from west to east, respectively. The western Atlantic samples comprised the Caribbean, Gulf of Mexico and southern Atlantic (Brazil), treating Rocas Atoll and St Peter and Paul separately. The allocation of sequences in different Caribbean subdivisions is based on the main regional patterns of marine currents as shown in population genetic studies of very common Caribbean reef species (Andras et al., 2013, purple sea fan; Truelove et al., 2017, lobster; Labastida-Estrada et al., 2019, lionfish). The exact allocation of sequences to populations is shown in Table 1. We also performed an AMOVA test to test for differentiation between western and eastern Caribbean locations as previous studies have indicated the presence of such population subdivision in the Caribbean region (Baums et al., 2006; Galindo et al., 2006). Here, the sequences were allocated in the two regions if they were based east or west of the Mona Channel and east or west of the Dutch Antilles. All analyses were performed using the Kimura 2-Parameter (K2P) nucleotide model of substitution (Kimura, 1980). The statistical significance of Φ-statistics was assessed against the null hypothesis of panmixia by 10,000 permutations of groups and haplotypes. Population differentiation was estimated using pairwise Φ_STs between populations (Weir & Cockerham, 1984) as implemented in Arlequin. We used the software MIGRATE v.4.4.4 (Beerli & Felsenstein, 1999; Beerli et al., 2019) to estimate the effective number of migrants (M_2→1, M_1→2). MIGRATE does not examine just the number of sequences but it takes into account in a Bayesian framework the relatedness of the sequences under the coalescence model. Because the dataset is limited and comprised of a single locus, we reduced the number of estimated parameters by comparing only the eastern (Mediterranean + Tenerife) vs western Atlantic (Wider Caribbean + Brazil) populations. The effective number of migrants per generation (N_em) was estimated between Western and Eastern Atlantic demes as θ•M. In the Bayesian search strategy available in MIGRATE, one long chain with four independent replicates was conducted for each run with 100,000,000 generations minus the first 10,000,000 steps which were discarded as burn-in. We used the following values for the heating chain: 1.0, 1.5, 3.0, and 1000,000.

Sequence divergence between populations was estimated in MEGA v.11 (Tamura et al., 2021) with the K2P distance method and p-distance for COI. Two methods were used to reveal the relationships among haplotypes: median-joining network analysis and the reconstruction of a phylogenetic tree. Genealogical analysis of the mitochondrial haplotypes was performed in PopART (Leigh & Bryant, 2015) with the median-joining network algorithm (Bandelt et al., 1999; ɛ = 0). Haplotype networks were built for all sampled regions together (see Table 1), for the Mediterranean locations (Ikaria, Crete and Malta separated) and one for the wider Caribbean (Caribbean and Gulf of Mexico separated). Phylogenetic analyses of COI haplotypes were conducted with the Maximum likelihood method as implemented in raxmlGUI 2.0 (Stamatakis, 2014; Edler et al., 2021). The program ModelTest-ng (Darriba et al., 2020) was used within raxmlGUI 2.0 to estimate the best model of nucleotide substitution for the COI haplotype alignment. Clade support was assessed with the transfer bootstrap expectation (Lemoine et al., 2018). The RaxML-ng algorithm (Kozlov et al., 2019) within raxmlGUI 2.0 was used with the following line of commands: raxml-ng –all –msa <FILENAME> --model TIM2 + G –prefix <FILENAME> --seed 954714 –bs-metric tbe –tree rand{10} –bs-trees 1000. The resulting tree file was uploaded to iTOL v6 (Letunic & Bork, 2021) for tree visualization.

Results

After cleaning and end-trimming, the final sequence alignments used for analysis consisted of 196 specimens for COI (502 bp) and 198 specimens for Cytb (338 bp). The final datasets consisted of new sequences and sequences downloaded from GenBank (N = 390 for COI; N = 198 for Cytb), as explained in the Materials and methods section. We identified 195 COI haplotypes, of which 147 were singletons and 118 haplotypes for Cytb, of which 94 were singletons. The highest number of haplotypes was observed in the Wider Caribbean and Caribbean in COI, 120 and 103, respectively (Table 2, Figure 2). The lowest values were observed in St Peter and Paul Archipelago (Table 2, Figure 2). The highest numbers of haplotypes (N = 72) for Cytb were identified in the Caribbean and the lowest (N = 20) in Ikaria of Greece (Table 3, Figure 3). The COI Haplotype 40 was the most common one (N = 68) and was exclusively found in Tenerife (Canary Islands, Spain) (N = 6) and in the Mediterranean (N = 62) (Figure 2). The COI Haplotype 1 was the second most common one (N = 18) and was found in the western Atlantic; in both the Caribbean (N = 14) and Brazil (N = 3), however, one specimen originated from the Mediterranean (Figure 2). The Cytb Haplotype 74 was most common in Tenerife (N = 9) and the Mediterranean (N = 10) (Figure 3). The second most common Cytb Haplotype 18 was found in the Caribbean (N = 16), but not in Tenerife or the Mediterranean Sea.

Fig. 2. Median-joining network based on COI haplotypes of Hermodice carunculata from the Caribbean, Brazil, Tenerife and the Mediterranean. The haplotypes are based on sequences generated in this study (GenBank accession # MF001520–- MF001715), as well as those from Ahrens et al. (2013; GenBank accession # KC017475-KC017594) and from Chatzigeorgiou et al. (2014; GenBank accession # KF878397-KF878476). The network was estimated with the median-joining network algorithm in PoPArt with epsilon = 0. Colours indicate the geographic origin of sequences, each circle is a different haplotype and the size of circles is proportional to the number of fireworms sharing this haplotype.

Fig. 3. Median-joining network based on Cytb haplotypes of Hermodice carunculata from the Caribbean, Tenerife and Mediterranean. The haplotypes are based on sequences generated in this study (GenBank accession # MF001716–- MF001913). The network was estimated with the median-joining network algorithm in PoPArt with epsilon = 0. Colours indicate the geographic origin of sequences, each circle is a different haplotype and the size of circles is proportional to the number of fireworms sharing this haplotype.

Table 2. Genetic diversity and summary statistics of Hermodice carunculata based on COI sequences

Geographic groupings of sampling sites are listed on Table 1. N, number of samples; H, number of haplotypes; n _e, effective number of haplotypes; S, segregation sites; π, nucleotide diversity; θ, Theta, Tajima's D and Fu's Fs. Standard deviation values in parentheses. Significant values are represented by ***, ** and * for P < 0.001, P < 0.01 and P < 0.05 (only for Tajima's D), respectively. Statistical significance for Fu's Fs: Not significant, P > 0.02.

Table 3. Genetic diversity and summary statistics of Hermodice carunculata based on Cytb sequences

N, number of samples; H, number of haplotypes; S, segregation sites; π, nucleotide diversity; θ, Theta; Tajima's D and Fu's Fs. Standard deviation values in parentheses. Significant values are represented by ***, ** and * for P < 0.001, P < 0.01 and P < 0.05, respectively. We used the Kimura 2-Parameter nucleotide model of substitution.

For the COI gene the Wider Caribbean was the region with the highest number of haplotypes both by simple count and by the estimation of the effective number of haplotypes (Table 2). The lowest effective number of haplotypes was observed in the Mediterranean despite the high number of sequenced fireworms (N = 138). The nucleotide diversity index π was highest (0.111) in Brazil and lowest (0.007) in Tenerife. In contrast, θ was highest in the Wider Caribbean (0.035) and lowest in the Mediterranean (0.008) (Table 2). The neutrality test Tajima's D was significantly negative in all locations, except for Brazil, indicating that populations may have undergone population expansions in the past. Fu's Fs statistic is a more sensitive indicator of population expansion and was significantly negative in all four data partitions. Strong deviations from neutrality equilibrium were expected because of the large number of singletons in every location (see Figures 2 and 3). For the Cytb gene, the nucleotide diversity index π was highest (0.015) in the Caribbean and lowest (0.009) in Tenerife. Furthermore, θ was highest in the Caribbean (0.0340) and lowest in Tenerife (0.025) (Table 3). Similar to the COI results, the neutrality tests Tajima's D and Fu's Fs were negative in all sampling locations, indicating population expansions of the fireworms in the past.

The COI haplotype network (Figure 2) indicates the presence of two distinct clusters of sequences, those from the Caribbean/Brazil and those from Tenerife and the eastern Mediterranean. Three of our sequences from Ikaria, Greece (MDGRHC23/Hap100, MDGRHC285/Hap170 and MDGRHC33/Hap172) and three published sequences (Ahrens et al., 2013) from Crete, Greece (MDGRCRKCOI17526/Hap1, MDGRCRKCOI17540/Hap45) and one from Malta (MDMTKCOI17553/Hap54) clustered with the Caribbean specimens (and the Brazilian specimens as seen in Figure 2). The region-specific haplotype networks (Supplementary Material 1 for the Mediterranean and 2 for the Caribbean + Gulf of Mexico) did not indicate any clear small-scale geographic patterns, justifying the pooling of individual worms from these regions in statistical tests. In the Mediterranean, 62 out of 138 H. carunculata from Crete (N = 42), Ikaria (N = 18) and Malta (N = 2) shared Haplotype 40 (Supplementary Material 1). Twenty-two of the haplotypes were singletons. In the wider Caribbean, the most common Haplotype (1) was present in 14 H. carunculata in the Caribbean and four in the Gulf of Mexico. The two areas shared four other haplotypes (16, 19, 29 and 66). One hundred singletons were recovered from both the Caribbean and the Gulf of Mexico (Supplementary Material 2).

Similarly, the Cytb haplotype network (Figure 3) shows two distinct clusters of haplotypes, the Caribbean one and the Tenerife/Ikaria-Greece one. A haplotype from Ikaria, Greece (GRHC33/Hap87) is clustering within the Caribbean haplotypes. The same sample is also clustering with the Caribbean haplotypes in the COI network. Two Caribbean samples from Guadeloupe (CRGPHC46/Hap 39) and one from Curaçao (CRCUHC178/Hap 8) clustered with the Tenerife/Ikaria-Greece haplotype group.

The Maximum likelihood analysis (Figure 4) yielded similar results as the haplotype network analysis. Almost all Hermodice sequences from the Mediterranean and Tenerife formed a well-supported clade (Figure 4, in blue; 85% bootstrap transfer values). Surprisingly, six haplotypes from the Mediterranean were placed among the Caribbean and Brazilian Hermodice (Figure 4, all six Mediterranean sequences denoted in bold blue). In fact, three Mediterranean sequences (Malta: MD_MT_KC017553 of Hap_54; Crete: MD_GR_CR_KC017526 of Hap 1; Ikaria: MD_GR_HC23 of Hap_100) were identical to sequences from the Caribbean and Brazilian specimens. The Caribbean or Brazilian Hermodice (Figure 4, in red) did not form distinct clades in the ML topology.

Fig. 4. Maximum likelihood consensus tree based on 195 COI haplotype sequences of Hermodice carunculata. Phylogenetic analysis was conducted with RAxML-ng as implemented in raxmlGUI 2.0. Tree was re-rooted at mid-point in iTol v6. Branches in red indicate haplotypes from the Caribbean and Brazil whereas branches in blue indicate haplotypes from the Mediterranean and Tenerife. Haplotypes 1, 45, 54, 100, 170 and 172 (in bold blue) include specimens from Mediterranean and Tenerife but are clustered with Caribbean and Brazilian specimens. Numbers above branches indicate bootstrap transfer values (only >70% values are depicted for clarity). BR, Brazil; BR CAR, shared haplotype between Caribbean and Brazil.

The AMOVA test based on COI indicated that there is significant population structure (Φ_ST = 0.644, P < 0.001), when partitioning the DNA sequences into four groups (Wider Caribbean, Brazil, Tenerife and the eastern Mediterranean) (Figure 1, Table 4). Most of the genetic variation (~62%) is allocated among groups (Table 4). The Φ_CT statistic was also highly significant (0.623, P ≪ 0.001), indicating that the within population component of the MOVA table was a significant source of variation (~35.6%, Table 4). The AMOVA test based on Cytb data indicated that there is significant population structure (Φ_ST = 0.584, P < 0.001), when partitioning the data in two groups (Caribbean, Tenerife + Ikaria-Greece) (Table 5) and in three groups (Caribbean, Tenerife, Ikaria-Greece; Φ_ST = 0.542, P ≪ 0.001) (data not shown). Pairwise Φ_ST comparisons among the three sampled locations (Caribbean, Tenerife and Ikaria-Greece) indicated that Caribbean H. carunculata are significantly distinct from both Tenerife (0.572) and Greece (0.563). Comparisons between Tenerife and Greece did not indicate significant differences (0.009). The AMOVA analysis (results not shown) between western and eastern Caribbean populations did not indicate the presence of significant population structure in H. carunculata (F _ST = 0.0036, P = 0.765).

Table 4. Analysis of molecular variance (AMOVA) of Hermodice carunculata based on COI sequences between the western and eastern Atlantic populations

The western Atlantic group was consisted of the Caribbean locations (Curacao, Panama, Florida/Bahamas, Puerto Rico, U.S. Virgin Islands, Guadeloupe, Dominica) including Flower Gardens, Sonnier Banks, Panama City, in Florida as well as Brazil (St Peter and Paul Archipelago and Rocas Atoll). The eastern Atlantic group was consisted of Tenerife and the Mediterranean locations Malta, Crete and Ikaria. Dominican Republic, Mexico and Rio de Janeiro were excluded from this analysis due to the small sample size.

F _ST = 0.644, P < 0.001; F _CT = 0.623, P ≪ 0.001; F _SC = 0.054, P ≪ 0.001.

Table 5. Analysis of molecular variance (AMOVA) of Hermodice carunculata based on Cytb sequences

Polychaetes were sampled from the Caribbean Sea, Tenerife (Canary Islands) and Ikaria (Greece).

F _ST = 0.584, P < 0.001; F _CT = 0.584, P = 0.33495; F _SC = 0.0006, P = 0.07475.

All pairwise Φ_ST comparisons among populations confirmed the presence of significant sequence divergence among most populations, especially those from the opposite sides of the Atlantic Ocean (COI; Table 6). All comparisons between the two sides of the Atlantic resulted in significant pairwise Φ_STs. The biggest differences occurred between eastern Atlantic and western Atlantic populations: Tenerife vs Panama City, Florida (0.749), and Tenerife vs St Peter and Paul Archipelago (0.735). The smallest pairwise Φ_ST values were observed among most of the comparisons of Caribbean populations, which were either 0 or very close to 0. The Brazilian population from St Peter and Paul Archipelago was different from all others, including those from Rocas Atoll. The polychaetes from Rocas Atoll were genetically much more similar to those of the Caribbean Sea than those from the St Peter and Paul Archipelago. The islands of Crete (eastern Mediterranean) and Tenerife (eastern Atlantic) displayed the biggest genetic difference (0.072), and the Mediterranean islands of Malta and Ikaria the smallest (0.024) (Table 6). Pairwise Φ_ST comparisons based on Cytb sequences among the three sampled locations (Caribbean, Tenerife and Ikaria of Greece) indicated that Caribbean H. carunculata are significantly distinct from both Tenerife (0.572) and Ikaria (0.563) (Table 7). Comparisons between Tenerife and Ikaria did not result in a significant Φ_ST value (0.009).

Table 6. A comparison of pairwise population differentiation of Hermodice carunculata based on COI sequences corrected by the Kimura 2-Parameter distance

Significant values (P < 0.05) in bold. Rio de Janeiro, Dominican Republic and Quintana Roo of Mexico were omitted from the pairwise comparisons because they are represented by ≤3 samples. The two samples from Statia were included in the samples from Guadeloupe. Peter Paul = St Peter and Paul Archipelago, Rocas At. = Rocas Atoll, FL_Bahamas = east Florida/Bahamas, Flower_Gar = Flower Gardens, Sonnier_B = Sonnier Banks, FL_PanamaC = Panama City, Florida.

Table 7. Average sequence divergence (%) between and within locations based on COI sequences as estimated in MEGA v. 11

Analyses were conducted using the Kimura 2-parameter model and in parentheses those with P-distance.

Using Migrate-N with the COI data, the mode of the estimated mutation-scaled effective population size of the western Atlantic is Θ = 0.0003 and of the eastern Atlantic Θ = 0.04410. The mutation-scaled size of the western Atlantic population is about 0.068% of eastern Atlantic H. carunculata. The 95% credibility interval for the effective population size of the western Atlantic population is 0–0.09993 and for the eastern Atlantic is 0.03107–0.05727. The estimated effective number of migrants per generation [N _em = (M_{i →j}* Θ_j)/4] between the group Wider Caribbean + Brazil and the group Med + Tenerife is 0.85 from west to east and 0.0000025 in the opposite direction.

The average sequence divergence based on K2P distance for the COI gene between the Caribbean specimens and those from the Mediterranean was 2.6, 2.8% compared with Tenerife and only 1.2% different from the Brazilian specimens (Table 7). Similar values were estimated when the sequence divergence analysis was conducted with P-distance. The pairwise nucleotide differences between populations for Cytb indicated that the Caribbean specimens were on average 2.1% divergent from fireworms collected from both Ikaria of Greece and Tenerife (Table 8).

Table 8. Average sequence divergence (%) between and within locations based on Cytb sequences as estimated in MEGA v. 11

Analyses were conducted using the Kimura 2-parameter model.

Discussion

The fireworm H. carunculata presents a very interesting case of a benthic omnivore with amphiatlantic distribution (Ahrens et al., 2013; Schulze et al., 2017). We investigated the genetic diversity, effective number of migrants and population structure of H. carunculata populations collected from three regions: the Caribbean, the temperate eastern Atlantic (Tenerife of Canary Islands) and the eastern Mediterranean Sea (Malta and the Greek islands of Crete and Ikaria). Since samples were collected over a wide geographic range, we expected to find significant sequence divergence in H. carunculata as previously reported by Ahrens et al. (2013) and Chatzigeorgiou et al. (2014). Indeed, data from two mitochondrial markers indicated that there is a 2–2.5% divergence between the amphiatlantic populations of H. carunculata, confirming our expectations. The presence of relatively low genetic divergence suggests that these allopatric populations belong to the same species and not to two species as suggested in Yáñez-Rivera & Salazar-Vallejo (2011). A step towards testing the species integrity of H. carunculata and the biological significance of the genetic differentiation among the amphiatlantic populations, would be to perform breeding experiments among western and eastern Atlantic specimens. If breeding experiments produce sexually reproductive offspring, the interbreeding specimens may belong to the same species, a criterion central to Mayr's (1942) biological species concept. Successful interbreeding, however, cannot exclude the possibility of the amphiatlantic fireworms being two true species that have not yet reached the stage of reproductive isolation (general lineage species concept; de Queiroz, 2007). Such experiments may be feasible since fireworms are hardy polychaetes (Ahrens et al., 2013) and may be transported live between labs for this purpose.

Hebert et al. (2003) estimated the average COI genetic divergence based on p-distance between 128 pairs of congeneric species of annelids to be 15.7% ± 6.2%, a value much higher than the one we reported here for H. carunculata. Comparisons between polychaete species pairs exhibited much higher sequence divergence values than most other metazoans (Hebert et al., 2003). There was considerable variation around the average sequence divergence among the species pairs of all 11 phyla examined by Hebert et al. (2003), but the vast majority of the metazoan species pairs (>98%; N = 13,320) showed >2% divergence. The amphiatlantic H. carunculata specimens differ by only 2.5–2.7% with P-distance and 2.6–2.8% with K2P distance (2.1% for Cytb) (Tables 7 and 8), suggesting that the observed divergence is most likely attributed to disconnected populations rather than two different species. We are aware that extrapolations of sequence divergence values between different species, even within the same family, are not warranted. Hudson and Turelli (2003) have strongly cautioned against neutrality-based predictions regarding mtDNA variation and divergence, emphasizing the stochastic nature of lineage sorting between species. Ideally, multiple genes could be used from the independently evolving mitochondrial and nuclear genome, since serious flaws have been identified when using small fragments of DNA to estimate the number of species (Hickerson et al., 2006; Rubinoff et al., 2006; Dasmahapatra et al., 2010). A limitation of our study is that both genes are in the mitochondrial genome and are inherited as a single locus. Reciprocal monophyly of mitochondrial and nuclear genes would strengthen our assertion about the historical separation of the amphiatlantic populations of H. carunculata.

The Atlantic basin is a major biogeographic barrier, especially for species with intertidal and sublittoral distribution (Briggs, 1974). Two reciprocally monophyletic amphiatlantic clades were found in the brooder polychaete Hediste diversicolor, with an average of 6.6% sequence divergence separating Canadian and European specimens (Breton et al., 2003), supporting the hypothesis that they represent two distinct species. However, there are several examples of polychaetes, including fireworms with amphiatlantic distribution with or without exhibiting significant population structure. The fireworm Eurythoe complanata, once regarded as a cosmopolitan species, consists of three deeply divergent lineages, one of which (Clade 1) is distributed from the Caribbean Sea to southern Brazil and islands off the west coast of Africa (Barroso et al., 2009). There was significant population structure among the Caribbean, Brazilian and African island populations of Clade 1, and to a lesser degree between Caribbean and Brazilian populations. A significant contributor to the population structure of E. complanata in the western Atlantic was the Amazon River, restricting gene flow between populations from north and south (Barroso et al., 2009). Even though we did not identify any significant differences between the Caribbean and the Brazilian population of Rocas Atoll, our analysis identified the St Peter and Paul Archipelago population as a genetically distinct one (Table 7). The fact that there is limited or no population structure between the Caribbean and the Brazilian samples from Rocas Atoll indicates the capacity of H. carunculata for long dispersal, probably due to its pelagic larval stage. Toso et al. (2020) studied the early development stages of H. carunculata in the laboratory and based on the duration of these first stages, they suggested the potential for a long planktotrophic period (i.e. teleplanic larvae); long enough to explain the absence of population structure between the Caribbean and Rocas Atoll as well as the small amounts of genetic divergence across the Atlantic Ocean. The linear distance separating Rocas Atoll and the edge of the Brazilian continental shelf is ~180 km, therefore the presence of teleplanic larvae and potential stepping-stone populations between these areas should facilitate the dispersal of fireworms along the South American coast to reach the Caribbean islands or vice versa.

The linear distance between the islands of St Peter and Paul Archipelago and the edge of the continental shelf of Brazil is 5 times longer compared with Rocas Atoll and forms an effective barrier for the trochophore larvae of H. carunculata. Within the Caribbean and the Mediterranean we did not observe significant population structure but this may be more indicative of the ability of the chosen mitochondrial marker to detect population subdivisions at this scale than the absence of fine structure. Microsatellite and SNP markers are more suitable for fine resolution structure studies within the Caribbean (e.g. Baums et al., 2006; Andras et al., 2013; Beltran et al., 2017) and should be applied to H. carunculata in order to test our hypothesis.

The suggestion that a coastal, common species such as H. carunculata is distributed across the Atlantic Ocean (Ahrens et al., 2013) is unsurprising especially in light of the long-lasting early larval stages observed by Toso et al. (2020). There are at least 134 opisthobranch gastropod species with amphiatlantic distribution, and 26.2% of those include Caribbean species (Garcia & Bertsch, 2009). Nineteen out of 22 species of prosobranch gastropods with teleplanic veligers are exhibiting amphiatlantic distribution (Scheltema, 1989). Although some of these gastropods may represent complexes of cryptic species, the Atlantic basin does not represent a very effective barrier to gene flow, especially for species with long-lasting planktonic larvae.

The inclusion of specimens from Tenerife to the Mediterranean-specific clade of H. carunculata shows that the Strait of Gibraltar does not impose restrictions on the gene flow between the eastern Atlantic (Canary Islands) and Mediterranean specimens and provides clues about the origins of these fireworm populations. In between these areas there is the Atlantic coast of Morocco and West Sahara, situated about 100 km away from the easternmost Canary Island of Fuerteventura which could serve as a stepping-stone. The Atlantic coast of Morocco might be a more likely origin of the connection between Mediterranean and Canary Island populations but unfortunately, this area was not represented in the current study. Patarnello et al. (2007) reviewed the significance of the Atlantic–Mediterranean transition in marine species concluding that the Strait of Gibraltar is a permeable biogeographic barrier, restrictive for some species but not for others. Multiple factors such as historical separations, glacial movements, and species' life histories have shaped the distribution patterns of fauna across the Strait of Gibraltar (Patarnello et al., 2007). The spill-over of eastern Atlantic fauna during the re-flooding of the desiccated Mediterranean basin about 5 Mya is a likely origin of the Mediterranean fireworms. Present day surface currents (Millot & Taupier-Letage, 2005) greatly favour the entering of organisms from the Atlantic Ocean to the Mediterranean Sea but impede the opposite movement. The geological history and the prevailing surface currents nowadays might contribute to explaining the greatly uneven genetic flow estimated in Migrate-N, with magnitudes of order higher estimated effective number of migrants per generation from the Atlantic Ocean to the Mediterranean Sea. The haplotype similarity between Tenerife and Mediterranean H. carunculata specimens indicates a common origin and genetic data from the Atlantic coast of Morocco could shed some further light on the biogeography of the region.

In a study within the eastern Mediterranean, Chatzigeorgiou et al. (2014) recovered two well-supported clades of H. carunculata: a Mediterranean-specific clade made of sequences from Crete and a clade consisting of both Cretan and Caribbean sequences. The Mediterranean-specific clade of Chatzigeorgiou et al. (2014) represents Clade I of Ahrens et al. (2013) and is similar to our clade, which contains mostly sequences from Greece, Malta and Tenerife. The same clade was recovered also by Righi et al. (2019) which included specimens from Porto Cesareo, eastern Italy. Clade I is not specific to the Mediterranean Sea any more as it contains sequences from Tenerife. The DNA sequences of this clade likely depict colonization events from eastern Atlantic sites; therefore neither the Almeria-Oran front nor the Siculo-Tunisian strait are effective biogeographic barriers for Hermodice carunculata as has been shown to be the case with other benthic invertebrates and fish (e.g. Rίos et al., 2002; Zardoya et al., 2004). The second clade recovered in both Chatzigeorgiou et al. (2014) and Righi et al. (2019) resembles our predominantly Caribbean group of sequences that includes very few sequences from the Mediterranean.

The presence of distinct clades within the same collection sites (Elounda and Alykes, northern Crete) in Chatzigeorgiou et al. (2014) is puzzling and suggests either the presence of naturally occurring sympatric but divergent lineages of H. carunculata or recent invasion events of ‘Caribbean’ specimens in the eastern Mediterranean. Similar patterns were observed by Ahrens et al. (2013), where in Clade Iia, a primarily Caribbean group included a haplotype mainly found in the Mediterranean, and a more geographically wide clade (Iic) contained haplotypes that were commonly found in Malta and Crete. Our analysis with new sequences supports the findings of Ahrens et al. (2013) and Chatzigeorgiou et al. (2014), as we have recovered several sequences of Mediterranean origin in the Caribbean clade and vice versa. Sympatric yet divergent lineages are not unusual and could be due to modern invasion events; e.g. Atlantic lineages of the barnacle Chthamalus proteus were found in Hawai'i (Zardus & Hadfield, 2005). Similarly, we suggest that these H. carunculata sequences provide possible evidence for anthropogenically induced dispersal. Arguably, the gene flow estimates between the western and eastern Atlantic populations as set up in the Migrate-N runs are driven by those ‘out-of-place’ specimens and they would have been even smaller if these specimens had been removed from the Migrate-N analysis. The highly asymmetric gene flow could also be an artefact of pulling differentiated populations together or the limited single locus data we are using. One locus may have limited information but Bayesian inference in Migrate-N may still achieve reliable results (Beerli, 2006). Polychaetes, adults and larvae, have been found in the ballast waters of ships (Chu et al., 1997); this dispersal mechanism is well known and has already been mentioned by Ahrens et al. (2013). Additional means of dispersal could be rafting on floating plastic (Barnes, 2002) or floating natural debris (Thiel & Gutow, 2005). The marine aquarium trade is another possible example of anthropogenically induced dispersal, where tons of ‘live rock’ substrata are transported from continent to continent every year. Fireworms are cryptic in nature and it is quite possible that specimens could be accidentally transported and released through the aquarium trade. Climate change-driven warming seas are also becoming an important dispersal factor of H. carunculata as Mediterranean populations are expanding northward along the Tyrrhenian and Adriatic Seas (Righi et al., 2020). The warming Mediterranean waters and the recent expansion of the Suez Canal (Fishelson, 2001; Galil et al., 2015) have attracted invasive tropical marine species from the Red Sea, which have established populations in the eastern Mediterranean; these species have been referred to as Lessepsian migrants. Fireworm specimens from the Red Sea could improve our understanding of the phylogeographic history of H. carunculata (Schulze et al., 2017) and test for the presence of Lessepsian migrants at the easternmost Mediterranean coasts.

Alternatively, the presence of distinct sympatric lineages within the same location may indicate the disruption of connectivity by ecological differentiation (Doebeli & Dieckmann, 2003; Prada & Hellberg, 2013). Divergent selection may favour seemingly sympatric lineages, if they live in contrasting ends of an environmental gradient, such as depth. In the marine environment, depth is one of the main drivers of ecological specializations and has been shown to limit gene flow, even at short distances of a few dozen metres (e.g. Carlon & Budd, 2002; Prada et al., 2008). Unfortunately, we do not know the precise habitat preference of fireworms since they are detectable only while feeding, outside of their cryptic habitat. More detailed ecological studies on fireworms may identify environmental gradients that distinguish the habitats the two sympatric lineages occupy and also provide improved density estimates of the populations. The mutation-scaled population size of the western Atlantic H. carunculata is much smaller than the eastern Atlantic H. carunculata as estimated in Migrate-N. Despite the very large credibility values of this estimate, it is worth investigating further with more detailed ecological studies and more genetic data since H. carunculata is an important omnivorous benthic invertebrate in its distribution range.

Surprisingly, our study and the other two similar ones on this topic (Ahrens et al., 2013; Chatzigeorgiou et al., 2014) recovered a few fireworm samples whose geographic origin did not match their genetic affinities, indicating reciprocal recent invasion events or divergent selection. We suggest that the discrimination between these two hypotheses could be another area of a wider research interest, regarding H. carunculata. We wish to extend this work and include a geographically wider sample (e.g. Cape Verde Islands, Azores, Atlantic coast of Morocco, Red Sea), where these additional populations may prove to be stepping stones or intermediate populations and may improve our comprehension of the history of this ecologically important species.