Multigene Molecular Systematics Confirm Species Status of Morphologically Convergent Pagurus Hermit Crabs

Introduction In spite of contemporary morphological taxonomy appraisals, apparent high morphological similarity raises uncertainty about the species status of certain Pagurus hermit crabs. This is exemplified between two European species, Pagurus excavatus (Herbst, 1791) and Pagurus alatus (Fabricius 1775), whose species status is still difficult to resolve using morphological criteria alone. Methodology/Principal Findings To address such ambiguities, we used combinations of Maximum Likelihood (ML) and Bayesian Inference (BI) methods to delineate species boundaries of P. alatus and P. excavatus and formulate an intermediate Pagurus phylogenetic hypothesis, based upon single and concatenated mitochondrial (cytochrome oxidase I [COI]) and nuclear (16S and 28s ribosomal RNA) gene partitions. The molecular data supported the species status of P. excavatus and P. alatus and also clearly resolved two divergent clades within hermit crabs from the Northeast Atlantic Ocean and the Mediterranean Sea. Conclusions/Significance Despite the abundance and prominent ecological role of hermit crabs, Pagurus, in North East Atlantic Ocean and Mediterranean Sea ecosystems, many important aspects of their taxonomy, biology, systematics and evolution remain poorly explored. The topologies presented here should be regarded as hypotheses that can be incorporated into the robust and integrated understanding of the systematic relationships within and between species of the genus Pagurus inhabiting the Northeast Atlantic Ocean and the Mediterranean Sea.


Introduction
Although hermit crabs are one of the most morphologically and ecologically diverse group of decapod crustaceans, their evolutionary history at many taxonomic levels is far from being resolved [1][2][3][4]. More than 1000 species, 127 genera and 6 families are currently reported for the superfamily Paguroidea [5], that inhabit diverse biotopes from intertidal to deep seas [6]. However, the true taxonomic and ecological heterogeneity associated with hermit crabs is likely to be underestimated because many species and lifehistories appear to be undescribed [4]. One of the most diverse groups belong to the family Paguridae, with species distributed widely through all oceans [5].
The genus Pagurus Fabricius, 1775 is considered to be ancient, with fossils being assigned to the genus as early as the lower Cretaceous [7] and from the Cenozoic [8]. Pagurus is the least complex of all paguroids, sharing a reduced or virtually nonexistent rostrum, no sexually modified appendages other than the regular female egg-bearing pleopods, and having no penis (or sexual tubes) [3]. Despite its comparatively morphological conservatism, Pagurus exhibits a high degree of species proliferation. Recently, 172 species (of which five are fossil records) are recognised [9] , that possess specialized adaptations for housing stability, relying upon gastropod shells for protection. Such commensalism has constrained morphological evolution over 150 million years [7,8] by requiring a decalcified asymmetrical abdomen capable of looping into gastropod shells. Despite the abundance of hermit crabs in the North East Atlantic and Mediterranean Sea [6] many important aspects of their taxonomy [10], biology [11][12][13][14][15][16][17][18], ecology [17,[19][20][21][22], systematics and evolution [3,8,23] are poorly documented.
Systematic problems remain among the Paguridae, such as the polyphyletic genus Pagurus Fabricius, 1775 [4]. To date, most systematic studies on hermit crabs have been based on morphology with relatively few studies utilising molecular tools to resolve species status [4,[24][25][26] or to determine phylogenetic relationships among major taxa [2,4,23]. Within the region of the Northeast Atlantic Ocean and Mediterranean Sea Pagurus is represented by 23 species [6] and high morphological similarity among some species has resulted in the recognition of two groups characteristic of the North Atlantic and Mediterranean Sea [10]. In some cases, morphology suggests very close relationships between congeners, raising uncertainty about their independent species status. Such a situation exists between the ''alatus'' group, Pagurus excavatus (Herbst, 1791) and Pagurus alatus (Fabricius, 1775) that is still difficult to resolve using morphology alone (Table 1) [6,10]. In spite of the prevailing taxonomic challenges, Ingle [10] has recognized both species, based mainly from the differences observed in the dorsal aspects of shield, antennular peduncle, the shape of the right cheliped, the length of larger pereiopod and male pleopods. Furthermore, due to a lack of life history studies there has always been considerable confusion regarding ongoing synonymies that are assigned among P. alatus, P. excavatus and P. variabilis (A. Milne-Edwards & Bouvier, 1892) [10,27]. As currently recognised, P. excavatus is distributed southwards from the southern part of the Bay of Biscay and into the Mediterranean Sea. Pagurus alatus extends northwards into Iceland waters, but remains sympatric with P. excavatus in many southern regions of the area [10]. Here, we use molecular phylogenetic analyses of mitochondrial cytochrome oxidase I (COI), mitochondrial 16S ribosomal RNA (16S), and nuclear 28S ribosomal RNA (28S) DNA partitions to reconstruct the systematics of selected Pagurus species, in order to make inferences on taxonomic status of P. alastus and P. excavatus.

Pagurus diversity
The variation in COI diversity was examined among 11 species (Table 2 and 3). For each species, one to six representative individuals were analysed (Table 4), and where possible, from different geographical areas, yielding a total of 46 sequences. No insertions, deletions, stop codons or sequences indicative of pseudogenes [28,29] were observed, and BLAST searches confirmed that the sequences corresponded to decapod mtDNA COI. There was also evidence for base composition bias in the sequences, notably a pronounced underrepresentation of guanine at the third codon positions (35.3% T; 18.1% C; 28.5% A; 18.1%G) a phenomenon commonly observed in metazoan mitochondrial [30]. The COI alignment contained a total of 513 bp with 177 variable characters, of which 170 were parsimony informative (33.13%). The high observed percentage of parsimony-informative character suggests that COI is sufficiently diverse for intrageneric phylogeny and clearly resolved all eleven Pagurus species examined in the present study ( Figure 1) [31,32]. The 11 species comprise a monophyletic clade ( Figure 1) with an average between genetic species distance of 17.01% (Table 2). Among the Pagurus species, P. acadianus exhibits the least genetic divergence to P. bernhardus (6.80%) and in contrast, P. arcuatus and P. excavatus exhibited the highest genetic distances (23.10%) ( Table 3). P. pubescens exhibited the highest average distance values (Table 3) with a range of 11.50-21.70% (see also Figure 1).

Phylogenetic relationships
Overall, the phylogenetic algorithms (ML, BI) resulted in congruent topologies, delimiting the designated true Pagurus species (Figure 1). The resulting molecular phylogeny agrees in several respects with the current morphologically based classification of all species. All analyses support the basal placement of P. longicarpus and identify two main clades ( Figure 1). In contrast, the relationships among inner clades of Pagurus were poorly resolved. Clade I is represented by four the Northeast Atlantic Ocean and Mediteranean Sea species and clade II by six species of the North Atlantic Ocean (East and West coasts) and Bering Sea specimens.

Molecular systematic assignments
The three independent genes revealed concordant phylogenetic differences between P. alatus and P.excavatus. Pagurus alatus is substantially divergent from P. excavatus, with a mean divergence of 14.9% and 5.1% for COI and 16S sequences respectively ( Table 5). The 16S alignment was more conserved than the COI partition yielding 84/462 variable characters, of which 75/462 were parsimony informative. The 16S sequences are AT-rich (71.86%), indicating a moderate compositional bias. The pattern of nucleotide substitution was also biased in favour of transitions over transversions, yielding a ts:tv = 1.1 and for 28S a ts:tv = 3.2.
In the BI analysis, systematic positions of five species were not stable (Figure 2 A), but the trees underpinned by the 16S, 28S and the concatenated data partitions were broadly congruent with the COI hypothesis ( Figure 1). The combined molecular analysis, Morphological selected characters by Ingle (1985) Dorsal aspect of shield and associated appendages Pagurus alatus (Fabricius 1775) Pagurus excavatus (Herbst 1791) Segment 1 of antennular peduncle, outer distal margin with one, two or sometimes three spines.
Segment 1 of antennular peduncle, outer distal margin without or with small obtuse spines at the most.
Outer dorso-lateral process of antenna (of large circa 80 mm SL, specimens) reaching just beyond distal margin of segment 4 and acicle reaching well beyond distal extremity of cornea; breadth of cornea slightly exceeding 1\2 length of eye.
Outer dorso-lateral process of antenna (of large specimens circa 80 mm SL) not reaching to distal margin of segment 4 and acicle reaching only to extremity of cornea; breadth of cornea slightly less than 1\2 length of eye.
Outer (particularly upper) surface of right cheliped palm not strikingly concave.
Outer upper (and sometimes lower) surface of right cheliped palm strikingly concave.
Words in bold are highlighting the solely differences between species and ''SL'' is the abbreviation of shield length. doi:10.1371/journal.pone.0028233.t001 based on strong posterior probabilities values (100%; in Figure 2, CON), clearly reveals the two independent lineages here discussed.

Pagurus diversity
Systematics work to date in Pagurus has primarily been in the area of morphological taxonomy, with few hypotheses presented regarding species relationships. Generating such inferences has been difficult because of generally remarkable similarities in morphology among members of the genus. These similarities are perhaps surprising given the degree of the genetic divergence observed in this study between the most phenotypically closest species, P. alatus and P. excavatus (14.8%). A similar pattern has been observed in penaied shrimps [33], porcelanids crabs [34], diogenid hermit crabs [35] and has been attributed to stabilizing selection on morphological/ecological characters, or that they are on independent evolutionary trajectories. Likewise, strong stabilizing selection acting on morphological differentiation has been associated with life as a commensal of thallassinidean shrimps, accompanied by neutral molecular divergence [36].
The smallest mean intraspecific divergence values observed (Table 3) are possibly underestimated, because samples were obtained from a single locality. Global-scale phylogeography surveys of COI sequence diversity have estimated average intraspecific diversity values of less than 1% within crustaceans, whereas interspecific values typically are greater than 4% among congeneric species [37][38][39] and especially among decapods that can exhibit congeneric divergence values greater than 15% [40].
Elsewhere, five species of the genus Pagurus from Sea of Japan exhibited lower levels of genetic identity when compared with the genera Metapenaeus and Penaeus [25]. Here, the high genetic diversity observed among the pagurid species is in line with the observed morphological variability in informal morphological groups among adults and larvae described by Ingle (1985) and McLaughlin and Gore (1988), respectively.

Phylogenetic relationships
Ingle [10] delineated two main groups of Northeast Atlantic Ocean and Mediterranean Sea Pagurus based on three adults, and additional larval morphological characters [10] that are fully congruent with the current molecular systematic analysis. Furthermore, P. armatus, P. ochotensis and P. bernhardus that were assumed to be most related, based on morphology [10], share the most basal position in clade II. In addition, all species from clade I and II agree with two distinct larvae groups described by McLaughlin and Gore [41] based on the characteristics and species assigned. The inconsistent phylogenetic position of P. bernhardus, observed in Mantelatto et al. [4] was still unresolved here (Figure 1), represented by the weak support of the two focal nodes (,50% bootstrap support) within clade II. However, the observed COI molecular divergence does support a sister group relationship between P. acadianus and P. bernhardus, a taxonomic relationship derived also from morphological comparisons [7,42]. In summary, the phylogenetic patterns observed for Pagurus are consistent with morphological groups established by Ingle [10] and McLaughlin and Gore [41] but further analyses would be required to establish the precise cladistic relationships within the genus as a whole.    Table 4 for complement information).Two major clades have been roman number-coded, I and II: represent two groups defined previously by Ingel (1985) based on adult and larvae morphological characters. doi:10.1371/journal.pone.0028233.g001

Molecular systematic assignments
Concordance across molecular and morphological characters provides a reliable indicator of longstanding evolutionary independence and consequently provides an operational criterion for species recognition [43]. In the present study, morphological differences between P. alatus and P. excavatus were evident when two morphological variants were compared directly (Table 1), supporting molecular findings that justify separation at the species level. It is important to mention that a restricted collection from one biogeographical region of a species with a wide distribution cannot represent the complete range of morphological variation that is often found in decapods [4,27,44]. We also know that phenotypic plasticity in morphological traits (especially organisms with commensal behaviours) might be strongly influenced by environmental factors between different areas [45][46][47]. The precise cladistic relationships among the five species were not stable (Figure 2, A), most likely due to lower taxon sampling and also the inability of COI to accurately resolve deep nodes. However, the combined molecular analysis (Figure 2, CON) clearly reveals the evidence for long-standing evolutionary independence between P. excavatus and P. alatus.

Conclusion
Molecular data have not been used before to investigate systematic relationships among Pagurus of the Northeast Atlantic Ocean and Mediterranean Sea. Despite the perceived limitations regarding the use of morphological characters for inferring evolutionary relationships among commensal species, our molecular data support the morphological taxonomy. Our data may indicate the possible existence of two monophyletic groups (clade I and II), supporting previous assertions based on larval and adult morphological criteria. However, the current data confirm the complexity of the relationships within Pagurus, highlighting the absence of complete and integrated morphological descriptions for the diverse and heterogeneous members of the genus. Since the present taxonomic and geographic coverage is incomplete, the topologies presented here should be regarded as working hypotheses. Pagurus have diversified into a wide variety of marine habitats and exemplify classic commensal, anti-predator evolutionary traits. Thus, the group provides an excellent model for studying the interplay between speciation, neutral molecular divergence and potential stabilising selection on body form.

Sampling
Twenty nine hermit crabs of six species, Pagurus alatus (Fabricius, 1775), P. bernhardus (Linnaeus, 1758), P. cuanensis (Bell, 1845), P. excavatus (Herbest, 1791), P. prideaux (Leach, 1815), and P. pubescens (Krøyer, 1838), were collected from the Portugal (Costa Algarvia, Costa de Prata, and Azores), United Kingdom (North Wales), Norway (Bear Island Slide and Svalbard) and Italy (Sicily) between 2006 and 2008. Specimens were harvested as by-catch from rough ground bottom trawls from INRB-IPIMAR, Sicilian fisheries survey and collected by physical searches. Species were identified (Ingle, 1985) prior to the excision and preservation of muscle tissue in 95% ethanol (stored at 220uC) and whole body storage at 220uC. In order to accurately assign specimens to currently accepted (female) species of P. alatus and P. excavatus we used morphological criteria based upon four main groups of characters: the dorsal aspects of shield (Table 1), antennular peduncle, the shape of the right cheliped, length of larger pereiopod and male pleopods (Ingle 1985).

DNA isolation, amplification and sequencing
Total genomic DNA was extracted from approximately 1 mm 3 of muscle tissue or whole legs for small specimens using the Chelex dry release method [48]. The COI gene was amplified with alternative sets of primers depending on PCR reaction success following by the protocol develop by Costa et al. [34]. All PCRs were performed in a 25 ml final volume containing 16 PCR buffer, 3-4 mM MgCl 2 , 0.1-0.3 mM dNTP, 1 U Taq polymerase, 5-10 pmol of each primer, and 2-10 ng of DNA template ( Table 6). The thermal cycling conditions are listed in Table 6 for each set of primers. Following amplification, PCR reactions were cleaned by incubation with 10 U Exonuclease I (New England Biolabs H) and 1 U Shrimp Alkaline Phosphate (Promega H) at 37uC for 1 hour, followed by heating at 80uC for 5 min. Samples were sequenced by Macrogen Inc. (South Korea) using an Applied BiosystemsH 3730 sequencer.
All empirically derived sequences were manually checked for ambiguities in CodonCode Aligner version 1.3.0, aligned using the ClustalX plugin embedded within Mega 4 [49], prior to manual quality control and megablast annotation.

Phylogenetic relationships
The most popular barcode marker COI is generally used to study close to moderately deep interspecific taxon relationships of crustaceans [50][51][52][53][54]. To provide a comprehensive sister-species coverage and assessment of interspecific variation, Pagurus COI sequences from GenBank were merged with our data ( Table 4). In that propose we use using Kimura 2-parameters (K2P) genetic distances within and among species implemented in Mega 4.1, and compared to literature data. Amino acid translations of the target genes were examined to ensure that no gaps or stop codons were present in the alignment. To identify phylogenetic groups among the resulting eleven putative Pagurus species, the 46 COI sequences comprising 513 bp was analyzed using Maximum Likelihood (ML) and Bayesian Inference (BI) phylogenetic reconstruction methods. The crab Macropodia longipes (Milne-Edwards & Bouvier, 1899) (Brachyura:Inachidae) and the most phylogenetically closest related species, Dardanus arrosor (Herbst, 1796) (Anomura: Diogenidae) were used as outgroups.
Since Bayesian posterior probability support values (bpp) can often be inflated for certain clades, relative to ML bootstrap values [55], we constructed trees using both Bayesian approaches and ML. For these searches, we set the substitution model parameters calculated by jModeltest in RAxML 7.0.4 [56] and in MrBayes 3.1 [57], respectively. Ten independent ML analyses were conducted using GTR+I+G with invariant sites (I) and gamma distributed rates (G) [58] (see below) as the model (using the GTR+CAT setting) with 4 categories of rate variation (500 bootstrap replicates were undertaken for estimation of node support) for each partition on combined data. In order to find the ML tree, 10 independent runs of RAxML 7.0.4 were conducted. ModelTest [59] identified the HKY+I+G model [60] as best indicated by Akaike Information Criterion (AIC) [61], however, since this model is not implemented in the current version of RAxML, the GTR+I+G model was selected as the closest matching alternative. In MrBayes, two independent Markov chain Monte Carlo (MCMC) analyses were run using four chains for 5610 6 generations with the initial 1 million generations (20%) cycles discarded as burn-in. To check that stationarity had been reached, we monitored the fluctuating value of the likelihood graphically with Tracer v1.4 [62]. Once the parameters reach stationarity, a 50% majority rule consensus tree was obtained from the remaining saved trees. The consensus tree was selected from the posterior distribution and visualized using FigTree V.1.0 (http:// tree.bio.ed.ac.uk/software/figtree/). Since substitution rates  among the four nucleotides and among different nucleotide sites in mitochondrial protein-coding genes as been reported [63][64][65], codon models of substitution TrNef+G [66], F81 [67] and HKY+G [60] were implemented in our BI analyses for first, second and third codon positions respectively.

Molecular systematic analyses
Use of nuclear genes in addition to mitochondrial genes adds to the number of independent markers in a dataset, thus increasing the chances to understand the systematic relationships between and within P. alatus and P. excavatus (Table 7). Here we analysed partial sequences of nuclear 28S (385 bp), mitochondrial 16S (462 bp), and the barcode region of COI (540 bp). The three gene regions were partitioned separately according to the previously determined model parameters (Table 8) in subsequent BI analyses. Gaps in 16S and 28S sequences were treated as a fifth characterstate. BI analysis was conducted for each gene data set and the concatenated partition (CON) with the three gene regions partitioned separately according to the previously determined model parameters (Table 8) as described before. To evaluate the range of intrageneric sequence identity found among Pagurus species, we compared pairwise distances of COI and 16S (Table 5).

Acknowledgments
Thanks are due to many colleagues and their respective institutions (Dr. Markos Aleksandro, Axel Barlow, Dr Nia Whiteley, Dr. Simon Webster, and Wendy Grail from Bangor University; Dr. Ana Cristina Costa and João Brum from Azores University; Dr. Debbie Bailie from Queens University; Dr. Marco Arculeo from University of Palermo; and Dr. Maria Włodarska-Kowalczuk from Institute of Oceanology PAS from Sopot) that collaborated through this study by making available some essential specimens, literature and sharing interesting scientific discussions. We thank Dr. Judite Alves and Alexandra Cartaxana for archiving and being responsible for the Crustacean collection in the Natural and History Museum of Portugal.