Highly contrasted population genetic structures in a host–parasite pair in the Caribbean Sea

Abstract Evolution and population genetic structure of marine species across the Caribbean Sea are shaped by two complex factors: the geological history and the present pattern of marine currents. Characterizing and comparing the genetic structures of codistributed species, such as host–parasite associations, allow discriminating the relative importance of environmental factors and life history traits that influenced gene flow and demographic events. Using microsatellite and Cytochrome Oxidase I markers, we investigated if a host–parasite pair (the heart urchin Meoma ventricosa and its parasitic pea crab Dissodactylus primitivus) exhibits comparable population genetic structures in the Caribbean Sea and how the observed patterns match connectivity regions from predictive models and other taxa. Highly contrasting patterns were found: the host showed genetic homogeneity across the whole studied area, whereas the parasite displayed significant differentiation at regional and local scales. The genetic diversity of the parasitic crabs (both in microsatellites and COI) was distributed in two main groups, Panama–Jamaica–St Croix on the one hand, and the South‐Eastern Caribbean on the other. At a smaller geographical scale, Panamanian and Jamaican parasite populations were genetically more similar, while more genetic differentiation was found within the Lesser Antilles. Both species showed a signature of population expansion during the Quaternary. Some results match predictive models or data from previous studies (e.g., the Western‐Eastern dichotomy in the parasite) while others do not (e.g., genetic differentiation within the Lesser Antilles). The sharp dissimilarity of genetic structure of these codistributed species outlines the importance of population expansion events and/or contrasted patterns of gene flow. This might be linked to differences in several life history traits such as fecundity (higher for the host), swimming capacity of larval stages (higher for the parasite), and habitat availability (higher for the host).


| INTRODUCTION
The history and dynamics of marine populations living in the Caribbean Sea have been shaped both by patterns of ocean circulation and geological events. Sea level fluctuations, related to the Quaternary climatic oscillations since 2.5 million years ago (Ma), changed the geography and ecology of the region. Eight climatic cycles have been recorded since 800,000 years ago (Pillans & Gibbard, 2012), notably the last glacial maximum (26 to 21 ka BP, Clark et al., 2009) with a sea level fall of ca. 150 m (Clark et al., 2009;Peltier, 2002;Peltier & Fairbanks, 2006). Such eustatic variations may have affected the distribution and the population genetic structure of extant organisms.
The present-day marine currents are characterized by three main systems namely the "Caribbean current," the "Antilles current," and a large eddy from Panama to Costa Rica (Lessios, Robertson, & Cubit, 1984;Gyory, Mariano, & Ryan, 2013; Figure 1). The speed and the direction of these currents (e.g., East to West along the Caribbean current, South to North along the Antilles current) may have implications for the genetic patterns among populations (e.g., direction of gene flow).
Integrating the Caribbean marine currents into an oceanographic model, four connectivity regions have been proposed (Cowen, Paris, & Srinivasan, 2006): Eastern Caribbean, Western Caribbean, Bahamas, and Panama-Colombia ( Figure 1). This regional pattern leads to the prediction of high dispersal potential of marine larvae within each region, but limited exchange across them . More recently, Kool, Paris, Andréfouët, and Cowen (2010) refined these connectivity regions, sometimes in weak agreement with geographic distances. The three new regions defined by this second model are the Lesser Antilles, Bahamas-Northern Cuba, and Panama-Nicaragua.
In both models, a break between western and eastern regions is predicted, and Jamaica is suggested as a stepping stone between them.
Comparing the genetic structures of codistributed species can disentangle the relative importance of common history, present-day ecology, and life history traits on their evolutionary history (e.g., Criscione, 2008;Kool et al., 2010). Parasite-host pairs are necessarily codistributed species with the distribution of parasites overlapping that of their specific habitat (the hosts). This is even more constrained if the set of host species is limited for a given parasite species (Poulin, 2007).
Characterizing and comparing the genetic structures of such codistributed species allow the identification of geographical barriers to dispersal (e.g., DeBiasse, Richards, Shivji, & Hellberg, 2016) or clarifies the contribution of landscape fragmentation to their phylogeographies (Rodelo-Urrego et al., 2013). Moreover, when interacting species have contrasted life histories, their comparison can also reveal which life history traits predominantly affect dispersal and population size among populations despite the shared environment (Criscione, 2008;Kochzius et al., 2009), and ultimately the co-evolutionary history of a given host-parasite association (Du Toit, Van Vuuren, Matthee, & Matthee, 2013). Finally, host-parasite costructure studies may help predict the potential for local adaptation by determining the relative dispersal rate between a host and its parasite (Greischar & Koskella, 2007).
Here, we aim to understand how populations of a marine hostparasite pair are genetically structured in the Caribbean. We studied the irregular sea urchin Meoma ventricosa and its parasitic pinnotherid crab Dissodactylus primitivus (De Bruyn, Rigaud, David, & De Ridder, 2009;Telford, 1982). Both species are endemic to the Caribbean Sea and to neighboring American coasts, from Florida down to Brazil (Alvarado, 2011;Chesher, 1969;Wirtz, de Melo, & De Grave, 2009).
Meoma ventricosa lives at depths of 1-200 m on soft substrates ranging from small coral pebbles to sandy or fine sediments (Chesher, 1969). Dissodactylus primitivus is an ectoparasite of M. ventricosa F I G U R E 1 Sampled sites across the Caribbean Sea with schematic pattern of the main currents. Labels 1 and 2, delimited by dashed frames, denote Panama-Nicaragua and Lesser Antilles "connectivity regions" according to Kool et al. (2010). Colors of the arrows denote differences in current speed (red: >20 cm/s, green: <20 cm/s) on which it reproduces, finds a shelter, and feeds (Telford, 1982).
Prevalence of parasitism is high, with 75%-100% of the sea urchins infected by 1-21 crabs (De Bruyn et al., 2009 and unpublished data).
The crab consumes host tegument and spines (Telford, 1982), which induce wounds on the sea urchin (De Bruyn et al., 2009 and D. primitivus have pelagic larvae and therefore are prone to dispersal by marine currents (Emlet, McEdward, & Strathmann, 1987;Pohle & Telford, 1983). However, the respective abundances and swimming abilities of the planktotrophic larvae of pinnotherid crabs and of sea urchins differ sharply: (1) Crab fecundity is thousands of times weaker than that of the sea urchin which might result in a lower dispersal and a lower population expansion capacity (Emlet et al., 1987;Jossart et al., 2014), (2) crab larvae are known to be better "swimmers" which should decrease drifting by marine currents (Metaxas, 2013;Yednock & Neigel, 2011), (3) habitat suitability (the sea urchin's body) is smaller for the parasite which should decrease recruitment rate. These differences in life history traits could cause incongruence in the population genetic patterns of these two partners.
We investigated the genetic variation of this host-parasite association to determine how past and recent ecological contexts shape the population structure of M. ventricosa and D. primitivus. Using partial sequences of Cytochrome Oxidase subunit I (COI) and microsatellites, we addressed the following questions: (1) Do host and parasite exhibit comparable genetic structures? (2) Do these structures correspond to the connectivity regions from predictive models in the Caribbean area? (3) What are the respective contributions of Quaternary sea levels fluctuations (Clark et al., 2009), of present pattern of marine currents, and of differences in life-history traits in explaining the crab and sea urchin demographic history or gene flow patterns?

| Collections
We sampled crabs and sea urchins between 2006 and 2013 at 14 sites (13 for sea urchins) (Table 1; Figure 1). These sites were situated at the Lesser Antilles (St Croix, Saint Barthélemy, Antigua, Guadeloupe, Martinique, Bequia, Canouan, Barbados), Greater Antilles (Jamaica), or Central America (Panama) (Figure 1). Samples were collected by SCUBA diving or snorkeling at depths ranging from 1 to 22 m (Table 1). Sea urchins were collected individually in plastic bags that were immediately tied up after collection.
Immediately after the dive, a sample of each sea urchin (3-4 spines) and all the crabs captured on each host were isolated, labeled, and preserved in pure ethanol.
The total numbers of specimens used for microsatellite analyzes were 327 sea urchins and 410 crabs (Table 1). For COI analyzes, we sequenced a total of 297 sea urchins and 309 crabs (Table 1).

| DNA extraction
We extracted DNA from one pereiopod of each crab using the Chelex resin method (see the detailed protocol in Jossart et al., 2014) and from two spines of each sea urchin using the Qiagen DNeasy Blood & Tissue Kit.
Using SAMOVA 2.0, we performed (for crabs) a Spatial Analysis of MOlecular VAriance (SAMOVA, Dupanloup, Schneider, & Excoffier, 2013). We calculated Φ CT for seven possible groupings (from 2 to 8) in order to find the grouping that maximizes the genetic variance among groups. Analysis of molecular variance (AMOVA, Excoffier, Smouse, & Quattro, 1992) was performed with Arlequin 3.5 (significance of Φ values was determined by a permutation test of 10,000 randomizations).
We tested isolation by distance (IBD) with a Mantel test (Φ ST vs. km), using the software Mantel 1.19 (life.bio.sunysb.edu/morph/softmult.html). The geographical distance corresponded to the shortest distance avoiding islands/strips of land and was calculated using the path tool in Google Earth. IBD was performed for the whole dataset and inside the "Lesser Antilles" connectivity region (excluding St Croix) defined by Kool et al. (2010) (see Figure 1).
We used three methods to verify the existence of population expansion. For these analyzes, locations showing no differentiation in other analyzes were pooled (see Section 3). First, using Arlequin 3.5, we calculated the following: (1) Fu's F S statistic, testing for an excess of recently evolved haplotypes in an expanding population compared with a stable population (Fu, 1997). The significance of the F S was determined by a permutation test using 10,000 randomizations.
(2) The sum of squared deviation (SSD) between the observed distribution of the number of nucleotide differences and the unimodal mismatch distribution expected from population expansion (Rogers, Fraley, Bamshad, Watkins, & Jorde, 1996;Schneider & Excoffier, 1999). SSD was also calculated to evaluate a potential spatial (range) expansion (Ray, Currat, & Excoffier, 2003). The significance of the observed mismatch was verified by a test of goodnessof-fit (10,000 bootstraps). For the spatial expansion analysis, we also estimated the time of expansion (Schenekar & Weiss, 2011).
(3) Past changes in effective population size were evaluated using the Extended Bayesian Skyline Plot approach (EBSP) in BEAST 2.4.4 (Bouckaert et al., 2014). For crabs, accurate estimates were not possible, because populations were differentiated and could not be pooled. For sea urchins, BEAST was performed for 2.10 7 iterations (10% of burnin), using a pairwise divergence rate of 1.52% per million years (Lessios, 2008) and HKY as the substitution model. Trace file was checked (including ESS values always > 200) using Tracer 1.6. The skyline plot was performed with an R script developed by the BEAST authors.
T A B L E 1 Sampling information including island/country, site, GPS coordinates, depth, year, and total number of samples for COI and microsatellite (SSR) analyzes

| Microsatellite data collection and analysis
For crabs, we used ten loci already used for studying Jamaican populations (Jossart et al., 2013(Jossart et al., , 2014. These loci were multiplex amplified and genotyped with an AB 3730 DNA Analyzer (see Jossart et al., 2013 for detailed protocol and primer sequences).
For M. ventricosa, we used eight microsatellite loci (Jossart, Geyer, & Lessios, 2015). Microsatellites were amplified in simplex according to the tagged primer-method and genotyped in an AB 3130XL Genetic Analyzer (see Jossart et al., 2015 for detailed protocol and primer sequences). We evaluated (using POWSIM 4.1, Ryman & Palm, 2006) that these microsatellites had a statistical power ( Using Genepop 4.2.2, the frequency of assumed null alleles was calculated, and linkage disequilibrium was tested for each locus pair within each species (Rousset, 2008). We used the software FSTAT 2.9.3.2 (Goudet, 1995) to estimate number of alleles and allelic richness (AR). Differences between sites in mean AR were tested using  Figure 1).
To infer the most probable number of genetic clusters (K), we used STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). For and allele frequencies correlated among populations. We determined the most likely value of K using the original method (described in the STRUCTURE manual) and the method of Evanno, Regnaut, and Goudet (2005) implemented in STRUCTURE HARVESTER (Earl & von-Holdt, 2012). Bar plots were created using the software DISTRUCT 1.1 (Rosenberg, 2004). In order to confirm STRUCTURE's assignments to genetic clusters in D. primitivus, we used BAPS 6.0 (Corander & Marttinen, 2006).
We also used the software divMigrate (Sundqvist, Keenan, Zackrisson, Prodöhl, & Kleinhans, 2016) to detect potential asymmetric gene flow between pairs of populations. For this analysis, the undifferentiated sites from the same island were pooled together.
Significance of asymmetry (10,000 bootstraps) was assed using the tool implemented in the software.

| COI data
In the crab D. primitivus, we found 39 haplotypes in 308 sequenced individuals. They were distributed in two divergent clades (18 haplotypes in clade A vs. 21 haplotypes in clade B) separated by 10 substitutions (Figure 2, Appendix S1 Bold values differ significantly from zero (based on 10,000 permutations). The alpha value was corrected by Benjamini-Yekutieli method (Narum, 2006) and were 0.0094 for D. primitivus and 0.0099 for M. ventricosa. See Table 1

| Microsatellite data
In both D. primitivus and M. ventricosa, no linkage disequilibrium was detected between pairs of loci in each population (600 and 364 pairwise comparisons, alpha was Benjamini-Yekutieli corrected to 0.0071 and 0.0076, respectively) (Narum, 2006 In M. ventricosa, STRUCTURE identified that the most probable number (K) of genetic clusters was one. In D. primitivus, the most probable K was four for the original method and two for the Evanno method.
The bar plot for K = 2 ( Figure 3) showed that most individuals from Panama Jamaica and St Croix were assigned to one genetic cluster. On the other hand, Martinique, Bequia, Canouan, and Barbados were highly associated with the other cluster, while those from remaining islands (St Barthélemy, Antigua, Guadeloupe) had more intermediate assignments. For K = 4 (Figure 3), the same situation was observed except that Barbados segregated in a single genetic cluster. The most probable K value in BAPS was 6, and the same subdivisions as those obtained with STRUCTURE were globally observed (Appendix S5).
divMigrate did not detect any asymmetric gene flow in M. ventricosa. In D. primitivus, several instances of asymmetric gene flow from Barbados to other islands (Guadeloupe, St Barthélémy, St Croix and Jamaica) were identified (Appendix S11).

| Contrast between host and parasite and its potential causes
The two interacting species exhibit highly contrasting genetic structures within the Caribbean Sea. The genetic diversity of the parasitic crab D. primitivus is structured between two main groups. One is mostly found in the western part of the Caribbean and the other one, in the eastern part. On the contrary, the sea urchin host M. ventricosa exhibits no genetic structure (either in mitochondrial and nuclear markers) across the entire investigated geographic area.
Both species exhibited signs of population expansion but a more recent expansion or a larger population size for the sea urchin host might explain the sharp contrast between the two species. It is likely that there are different magnitudes of gene flow taking into account the dissimilar dispersal abilities of the two species. Whereas both adult crabs and sea urchins are able to move, such movements are very local (tens of meters). Therefore, dispersal capacity is related to pelagic larval stages. It is conventional to consider pelagic larval duration (PLD) as the main contributor to dispersal distance, although the relationship of PLD with genetic structure varies between species (Dawson, 2014;Faurby & Barber, 2012;Shanks, 2009;Shulman & Bermingham, 1995). The PLD of M. ventricosa is unknown while the one of D. primitivus is approximately 2 weeks (Pohle & Telford, 1983). Based on other tropical sea urchins with pluteus larvae (Emlet et al., 1987), it is probable that the PLD of M. ventricosa is at least equal to the one of D. primitivus. PLD might be one of the contrasting life history traits between these species but this need to be evaluated for M. ventricosa. At least three other factors might be linked to the contrasting dispersal abilities of M. ventricosa and D. primitivus. First, high fecundity has a positive influence on dispersal by increasing the number of potential migrants (Johnston, Miller, & Baums, 2012;Palumbi, 1994).
Third, there are differences in the availability and suitability of settlement habitats (Cowen & Sponaugle, 2009;Treml, Halpin, Urban, & Pratson, 2008). The distribution of sandy banks-favorable habitats for the sea urchin-can cover several square kilometers and are present in most areas of the Caribbean, allowing frequent settlement after longdistance dispersal for Meoma larvae. The suitable area for recruitment is much more limited for crab larvae. They must find a habitat populated by the sea urchin, locate a host, and settle on its body (an area much smaller than the sand patches on which sea urchin larvae can settle).
If sharp differences in gene flow exist between these species, it would have implications for the evolution of parasitic interactions.
Theoretical models (e.g., Blanquart, Kaltz, Nuismer, & Gandon, 2013;Gandon, 2002) and a meta-analysis (Greischar & Koskella, 2007) suggest that when gene flow is higher in the host than in the parasite, natural selection will act on more alleles in the host which would adapt faster thus limiting the impact of its parasites. On the other hand, very high levels of gene flow may prevent local adaptation (Lenormand, 2002

| Comparisons with other taxa and with proposed biogeographic regions in the Caribbean
Genetic homogeneity of M. ventricosa populations was observed across the Panamanian region and the Eastern Caribbean. This pattern was previously detected within the Caribbean in other taxonomic groups (Johnston et al., 2012;Purcell et al., 2006;Silberman et al., 1994), as well as in other sea urchins (Lessios, Kane, & Robertson, 2003;Lessios, Kessing, & Pearse, 2001;Lessios et al., 2012;McCartney, Keller, & Lessios, 2000;Zigler & Lessios, 2004). The present study confirms this last observation, not only for mitochondrial DNA, but also for microsatellites.
The observed pattern of D. primitivus indicated that populations in Panama-Jamaica, but also to a lesser extent St. Croix, were differentiated from those in the other islands of the Lesser Antilles.

| Refinements in the genetic structure of the parasitic crab
The crab populations can be considered as genetically homogenous across several sets of locations. This is the case for all sites from the same islands, confirming previous results obtained for Jamaican coasts (Jossart et al., 2013). We also observed a genetic proximity  F I G U R E 3 STRUCTURE bar plots for K = 2 (top), K = 4 (bottom) in Dissodactylus primitivus. Each line corresponds to an individual that was assigned with a certain probability to each genetic cluster Armstrong, 1994;Luckenbach & Orth, 1992;Yednock & Neigel, 2011).
However, in the zone South to Guadeloupe, the Caribbean east-west current has higher speed, promoting high gene flow between these southern islands (Martinique, Grenadines).
Finally, we observed a moderate segregation of Barbados population (in both COI and microsatellite analyzes) despite its geographical proximity to other islands. Moreover, some cases of gene flow calculated between Barbados and its neighboring islands were asymmetric.
This corroborates the study of Roberts (1997), who estimated that the larval import in coral reefs of Barbados is one of the lowest in the Caribbean. However, as genetic diversity of D. primitivus at this island is not lower than in other islands, an external input should be consid-

| Quaternary climatic oscillations
Our results suggest that Quaternary climatic oscillations (glacialinterglacial periods) Peltier & Fairbanks, 2006) had an influence on the distribution of genetic diversity of crab and sea urchin populations. Several sea level falls of 100-150 m below the present level were related to the succession of glacial episodes (starting 2.7 Ma, more important 0.8 Ma when the glacial periods became much stronger, and ending with the last glacial maximum 0.02 Ma; Pillans & Gibbard, 2012). This temporal pattern coincides with both the expansions calculated from mismatch and EBSP analyzes. A similar population expansion was also suggested for two Caribbean squirrelfishes (Bowen, Bass, Muss, Carlin, & Robertson, 2006), for a corallivorous mollusk (Johnston et al., 2012), for a sea urchin (Lessios et al., 2003), and for a pea crab (Ocampo et al., 2013). This expansion might be linked to an increase in habitat availability during interglacial periods. Indeed, M. ventricosa is a coralassociated organism that could have thrived with the coral reef expansion during these interglacial times (Baums, Scott Godwin, Franklin, Carlon, & Toonen, 2013;Bowen et al., 2006;Johnston et al., 2012). For the crab, population expansions were also suggested by mismatch analyzes. However, the large confidence intervals on the time estimates do not clarify whether these expansions happened synchronously everywhere and simultaneously with the host expansion.