1 INTRODUCTION

Hermit crabs are anomurans decapod crustaceans of the superfamily Paguroidea, which included more than 1 100 species (de Grave et al., 2009; McLaughlin et al., 2010). Most of the hermit crabs have an asymmetrical abdomen secluding in a shell of a gastropod. Although hermit crabs have colonized a range of inhabits from terrestrial to the deep sea (Tsang et al., 2011), they are best known from intertidal areas. As active carrion scavengers, hermit crabs are ecologically important groups of zoobenthos. (Williams and McDermott, 2004). Hermit crabs also represent an intermediate group of crustaceans from Macrura to Brachyura, which occupys an important position in crustacean evolution (Sha et al., 2015). Diogenidae Ortmann, 1892 is the second largest family of the superfamily Paguroidea. Most species of Diogenidae live in the ocean, and only one species has been reported from freshwater so far (See McLaughlin and Murray, 1990). Clibanarius Dana, 1852 is one of the speciose genera of Diogenidae and usually inhabits intertidal zones. It currently contains 59 species and 43 of them are known from the Indo-West Pacific (McLaughlin et al., 2007, 2010; Malay et al., 2018). They are commonly found in tropical or subtropical shallow water inhabiting sandy, estuary, mangroves, seagrass beds, rocky shore, and coral reef (Hirose et al., 2010; Sha et al., 2015). Most species of Clibanarius are small and share a superficial resemblance in morphology, increasing the difficulty of identification. The color patterns in the life of most Clibanarius species is often used as a diagnostic character for the species identification (Sha et al., 2015). It has been reported, however, that color differences within samples may not represent the different species (Knowlton and Mills, 1992; Hirose et al., 2010). In addition, the body color of hermit crabs easily fades away in alcohol, which makes identification more difficult.

Taxonomic groupings of Clibanarius have been proposed by Asakura (2005) based on the relative length of dactyls and propodi. According to Asakura's proposal, the species of Clibanarius can be divided into two morphological groups: (ⅰ) species have the dactyls shorter than the propodi of the second and third pereopods; (ⅱ) species have the dactyls equal to or slightly longer than the propodi of the second and third pereopods (Asakura, 2005). However, this taxonomic division has not been verified based on molecular data. In contrast to the wealth information on the phylogenetic relationships among the major taxa of Anomura (e.g., Ahyong et al., 2009; Tsang et al., 2011; Bracken-Grissom et al., 2013), few studies have focused on the interspecific relationships within Clibanarius (Hirose et al., 2010; Negri et al., 2014; Yoshikawa et al., 2018). Hirose et al. (2010) performed a preliminary phylogenetic analysis of Clibanarius based on partial sequences of the mitochondrial cytochrome c oxidase subunit Ⅰ (COI) gene, in which poorly resolved phylogenies prevented phylogenetic relationships of Clibanarius from an in-depth investigation. The close association of C. humilis, C. merguiensis and C. englaucus and the close relationships of C. longitarsus and C. striolatus were strongly supported. Negri et al. (2014) undertook phylogenetic analyses of selected Clibanarius species using the barcode region of the COI gene to distinguish C. symmetricus from C. vittatus. However, most nodes of the phylogenetic tree showed low support values. They only showed the close relationships between C. symmetricus, C. vittatus, C. lineatus, and C. sclopetarius. Yoshikawa et al. (2018) identified three specimens as Clibanarius virescens that had unknown coloration for Clibanarius based on the mitochondrial COI and nuclear histone H3 genes, in which phylogenetic relationships of Clibanarius species had received little attention. Overall, the interspecific relationship of Clibanarius is remained poorly resolved and the hypothesis of Asakura (2005) remained no attention.

It should be noted that the above-mentioned molecular phylogenies are based on analysis of a limited amount of sequence data, mainly mitochondrial DNA sequences, which could not provide enough phylogenetic resolutions. Two nuclear protein-coding genes, sodium-potassium ATPase α-subunit (NaK) and phosphoenolpyruvate carboxykinase (PEPCK), play elementary roles in multitudinous life forms and maintain good conservation throughout evolution (Tsang et al., 2008). These two genes have been widely used to recovery the phylogenetic relationships among higher taxonomic category in insects for many years (Friedlander et al., 1996; Wiegmann et al., 2000; Leys et al., 2002; Danforth et al., 2004, 2011). Recently, these two nuclear protein-coding genes have been increasingly utilized to resolve the phylogeny of the caridean shrimps (Li et al., 2011) and Anomura (Tsang et al., 2011). The mitochondrial DNA reveals speedy substitution saturation that restricts their practicability in resolving deep nodes, whereas the nuclear proteincoding genes are easy to align and provide a rich source of information that promotes a more precise recovery of phylogeny (Tsang et al., 2008). Therefore, the combination of mitochondrial and nuclear genes which has successfully resolved some controversial phylogenetic relationships is expected to provide more accurate phylogenetic information (Kiwaidae, Roterman et al., 2018).

To investigated phylogenetic relationships of 11 Indo-West Pacific Clibanarius species, both mitochondrial (16S rRNA and COI) and nuclear (NaK and PEPCK) genes are sequenced. The present study aims to examine whether phylogenetic groupings based on molecular data support or refute previously proposed taxonomic divisions and to elucidate the interspecific relationships within Clibanarius. Furthermore, we also aim to compare the practicability of mitochondrial genes or nuclear protein-coding genes in exploring interspecific relationships of Clibanarius. Establishing phylogenetic relationships of Clibanarius species will improve the understanding of taxonomic classification and evolutionary relationships of Anomura.

2 METERIAL AND METHOD 2.1 Sample collection

A total of 11 representatives of Clibanarius Dana, 1852 are utilized in this study, (detailed sampling information shown in Table 1). Species of the specimens are identified morphologically according to McLaughlin et al. (2007) and Sha et al. (2015). Among them, Clibanarius corallinus (Milne Edwards, 1848), C. virescens (Krauss, 1843), C. merguiensis de Man, 1888, C. englaucus Ball and Haig, 1972 and C. rutilus Rahayu, 1999 have ambulatory dactyl shorter than the propodus, while C. longitarsus (de Haan, 1833-1850), C. striolatus Dana, 1852, C. infraspinatus (Hilgendorf, 1869), C. eurysternus (Hilgendorf, 1879), C. humilis (Dana, 1851), and C. snelliusi Buitendijk, 1937 have ambulatory dactyl equal to or slightly longer than the propodus. Specimens used in this study are preserved in the Marine Biological Museum (MBM) in the Institute of Oceanology, Chinese Academy of Sciences, Qingdao (IOCAS). The outgroup taxa, namely Pagurus ochotensis Brandt, 1851 and Hippa adactyla (Fabricius, 1787), are chosen to root the phylogenetic trees. All the samples are preserved in 75%–95% ethanol.

2.2 DNA extraction, PCR amplification, and sequencing

The pereopod or cheliped of each hermit crab (about 15 mg) is used for DNA extraction using a QIAGEN DNeasy Blood and Tissue Kit. The total genomic DNA is eluted in 200 μL of sterile distilled H2O and stored in -20℃ freezer. Partial fragments from four genes, 16S rRNA (~500 bp), COI (~660 bp), NaK (~590 bp), and PEPCK (~590 bp), are amplified by polymerase chain reaction (PCR). The amplifications are executed in 25 μL volumes that incorporated 1–2 μL of template DNA, 1 μL forward and reverse primers, PCR Mix containing 1 μmol/L RED Taq, 22 mmol/L Tris-HCL, 55 mmol/L KCl, 1.65 mmol/L MgCl₂, and 220 μmol/L dNTP. The primers used are revealed in Table 2. For 16S rRNA and COI gene segments, the program cycling conditions are as follows: initial denaturation for 3 min at 94℃, followed by 33 cycles of 30 s at 94℃, 30 s at 50℃, and 50 s at 72℃, with a final extension at 72℃ for 3 min. The following primer combinations are used to amplify NaK and PEPCK gene segments: NaK for-a/NaK rev3, NaK for-b/NaK rev3, NaK for-b/NaK rev and PEPCK for/PEPCK rev, respectively (Tsang et al., 2008, 2014). The PCR profile for NaK is as follows: 3 min at 94℃ for initial denaturation, then 33 cycles of 30 s at 94℃, 30 s at 50–55℃ depending on individual samples, 50 s at 72℃ with a final extension at 72℃ for 3 min. The same program cycling conditions are employed for PEPCK with the extended time of 1 min 30 s at 72℃. The PCR products are purified by using of QIAquick gel purification kit (Qiagen) before sequencing. Purified PCR products are bidirectionally sequenced with ABI 3730xl DNA Analyzer (Applied Biosystems, Foster City, CA, USA) using the same primer pairs for PCR amplification. The sequences of the four genes of Hippa adactyla are obtained from GenBank. Similarly, NaK sequences of C. virescens-2, C. humilis-1, C. englaucus-2 and PEPCK sequences of C. virescens-2 and C. humilis-2 are acquired from GenBank. All sequences generated in the present study have been accessioned in GenBank, and corresponding accession numbers (MH508045– MH508087, MK076123–MK076149) are given in Table 1.

2.3 Data analysis

The chromatograms of forward and reverse sequence fragments are checked and assembled using SeqMan in the DNASTAR Lasergene software package. The inferred amino acid sequences are obtained to detect the presence of nonsense mutation using the program MEGA 6.0 (Tamura et al., 2013). Molecular diversity indices, such as the number of polymorphic sites, parsimony information sites, the insertions or deletions of nucleotide sequences, are achieved by Arlequin 3.5 (Excoffier and Lischer, 2010). PAUP*4.0b10 (Swofford, 2002) is used to examine the base compositional homogeneity within the taxa implemented by the Chi-square test. The sequences of the four genes are aligned separately using MUSCLE default settings that integrated into MEGA 6.0. After alignment, we splice the nucleotide sequences of the four/two genes to form a combined dataset for downstream phylogenetic analyses.

The phylogenetic trees are reconstructed using Maximum likelihood (ML) and Bayesian inference (BI) algorithms. The ML analysis is accomplished on the ATGC bioinformatics platform with PhyML 3.0 (Guindon et al., 2010), and Bayesian inference (BI) is conducted using MrBayes 3.2 (Huelsenbeck and Ronquist, 2001). In total, we analyze three datasets: all four concatenated genes, two concatenated mitochondrial genes, and two concatenated nuclear genes. The best-fit models of nucleotide substitution of each gene are selected using jModelTest v2.1.6 (Darriba et al., 2012). The best-fit model of the four concatenated genes and the concatenated two mitochondrial genes is GTR+I+G. The optional model of the concatenated two nuclear genes is TIM2+I+G. The best-fit models shown by the jModelTest for the 16S rRNA, COI, and NaK, PEPCK are TIM3+I+G, GTR+I+G, TIM1+G, and TIM2ef+G, respectively.

We conduct the ML analysis using the best-fit models of DNA substitution for the combined mitochondrial genes, nuclear genes, and all the four genes. For each analysis, if the model selected by the jModelTest is not available in the MrBayes and PhyML 3.0, we would use the most similar model to instead. The supporting value for each clade is estimated by 1 000 bootstrap replicates (BP). In a hybrid model, the partitioning of sequence data that come from multiple genes can provide a better estimation of phylogenetic relationships (Castoe et al., 2004; Brandley et al., 2005; Brown and Lemmon, 2007). Thus, we undertake the partitioned Bayesian analyses for each concatenated dataset, allowing specific substitution models to be assigned to each gene. A metropolis coupled Monte Carlo Markov Chains (MCMCMC) sampling is performed with four chains running simultaneous for 10 000 000 generations with trees sampled every 1 000 generations. The average standard deviations of the split frequencies (< 0.01) and potential scale reduction factor (RSRF, close to 1.0) are used as the criterions to validate the convergence of the analysis. After discarding the first 25% of the generations as burn-in, the 50% majority-rule consensus trees from the remaining trees are constructed to estimate posterior probabilities (PP) for each clade.

3 RESULT 3.1 Sequence characteristics

The combined dataset of the four gene fragments consists of 2 304 nucleotide sites. The 16S rRNA and COI sequences show a strong AT bias, with mean overall values of 68.8% and 61.1%, respectively. Amongst the 519 sites of 16S rRNA, 200 are variable and 152 are parsimony informative. Amongst the 654 sites in COI, 298 characters are variable, 241 are parsimony informative. Sequences of PEPCK are slightly G/C biased (54.4%), whilst the NaK is also slightly G/C biased (57.0%). The aligned sequences of NaK gene includes 582 nucleotide positions with 133 variable sites and 63 parsimony informative sites. The PEPCK gene includes 549 positions in which 158 are variable and 75 are parsimony informative. No stop codons or indel are observed from the COI, NaK, and PEPCK. A total of 31-bp deletion is detected in the 16S rRNA sequence alignment. There is no significant base heterogeneity across all codon positions in the three protein-coding genes (Chisquare P>0.05), except in the third position of the COI gene (Table 3).

3.2 Phylogenetic inference

Unexpectedly, a discrepancy between the results of mitochondrial and nuclear DNA datasets is found in the present study (Fig. 1). Moreover, the two analytic methods (ML and BI) yield different results for each concatenated dataset. The difference between the tree topologies received by ML and BI methods of the mitochondrial concatenated dataset is the position of the species C. snelliusi and C. eurysternus. These two species cluster together and form a clade with C. striolatus, C. longitarsus, C. infraspinatus and C. rutilus using the ML method (BP=51%), while form a clade with C. corallinus, C. virescens, C. merguiensis, C. englaucus, and C. humilis using the BI method (PP < 0.50). The difference between the tree topologies received by ML and BI methods of the nuclear concatenated dataset is the position of C. englaucus-1. C. englaucus-1 and C. longitarsus occupy the basal position in the ML tree (BP < 50%), while C. englaucus-1 cluster with C. merguiensis-1, C. englaucus-2, and C. humilis in the BI tree (PP=0.82). Nevertheless, all the analyses reveal consistently the close relationship among C. merguiensis, C. englaucus, and C. humilis. On the other hand, the tree topologies of the four-gene concatenated dataset are congruent regardless of the methods analyzed. Here, only the BI tree is presented with support values for both ML and BI analyses (Fig. 2). The interspecific phylogenetic relationships of Clibanarius are better resolved by using all four genes concatenated dataset, which could not be clearly uncovered in the other two datasets with lower support values. For these reasons, we use this consensus tree as the phylogenetic tree of Clibanarius.

As shown in Fig. 2, the tree topologies derived from the ML and BI analyses are highly congruent and some deep nodes are strongly supported. The clustering of C. rutilus with the dactyl shorter than the propodus and C. infraspinatus with the dactyl longer than the propodus is strongly supported by the analyses (BP=100%, PP=1.00) and form a clade with C. snelliusi with high support values (BP=91%, PP=1.00). C. longitarsus and C. striolatus cluster together with relatively high support values (BP=69%, PP=0.99), while the phylogenetic trees show relatively low support values for their relationship with C. rutilus, C. infraspinatus and C. snelliusi (BP=52%, PP=0.90). The branching positions of C. eurysternus and C. corallinus are not obtained in the current tree topologies. C. virescens represents an early divergence from C. merguiensis, C. englaucus and C. humilis with high support values (BP=90%, PP=0.99). Therein, C. merguiensis and C. englaucus with the dactyl shorter than the propodus clustered as a sister taxa with high supports (BP=93%, PP=1.00), forming a clade with C. humilis with the dactyl equal or subequal to the propodus (BP=99%, PP=1.00). It is worth noting that the supports values of some nodes for the interspecific relationships in the molecular tree are relatively low and further studies with more taxa sampling and genetic information are required to resolve this issue.

4 DISCUSSION 4.1 Implication of comparative studies on Clibanarius

This is the first study in which we investigated the phylogenetic relationship of Clibanarius based on both mitochondrial and nuclear DNA sequences. Previous studies on Clibanarius or other genera of diogenid phylogeny were based on morphological characteristics or/and few molecular data that relied on a single mitochondrial gene (Mantelatto et al., 2006; Tirelli et al., 2008, 2010; Hirose et al., 2010; Negri et al., 2014). Therein, they all showed lower support values and discrepancies in the critical nodes of Clibanarius phylogeny. In this study, we have explored the potential impact of the amount of sequence as well as the inclusion of multiple genes with distinct evolutionary histories on the phylogenetic reconstructions of Clibanarius. The lower resolution of interspecific phylogenetic relationships of Clibanarius is observed when mitochondrial DNA and nuclear DNA dataset are analyzed separately. Therein, two analyses methods product different tree topologies by using same dataset. By contrast, most of the interspecific relationships are robust and coincident in both ML and BI analyses based on the combined mitochondrial and nuclear DNA dataset. It has been recommended that long DNA sequences can yield a sufficient number of characters for robust phylogenetic inference and provide high phylogenetic resolution of relationships among closely related congeneric species (Minegishi et al., 2005; Cheng et al., 2011). The relatively short DNA sequences used for phylogeny construction may provide an explanation for the inconsistent phylogenetic patterns observed in previous reconstructions of Clibanarius phylogeny (Hirose et al., 2010; Negri et al., 2014).

In this study, we also demonstrate the importance of including different gene regions of molecular markers in the phylogenetic resolution of relationships among Clibanarius species. The characteristics of high copy number, rapid evolution rate and maternally inherited of mitochondrial genes makes them widely used as a molecular marker in phylogenetic studies (Moritz et al., 1987; Toon et al., 2009). Nevertheless, the rapid rate of nucleotide substitution among mitochondrial genes has limited their utility in higher levels of divergence (Toon et al., 2009). On the contrary, some of the nuclear protein-coding genes have the advantage of lower nucleotide substitution rate and are easy to align, which exhibit better performance in resolving phylogenetic relationships at deeper taxonomic levels (Moriyama and Powell, 1997; Chu et al., 2009). Therefore, a multigene phylogenetic approach is needed to construct a phylogeny, with slower-evolving genes settling deeper relationships and faster-evolving genes settling closer relationships (Palero and Crandall, 2009). It has been reported that a well-resolved topology is obtained to reconstruct phylogenetic relationships of the infraorder Stenopodidea (Chen et al., 2016) using mitochondrial and nuclear DNA sequences from a combined perspective. We consider it more likely that different mechanisms of nucleotide substitution and different evolutionary properties may explain the discrepancy between mitochondrial and nuclear DNA phylogenies of Clibanarius and suggest the significance of considering different gene regions of molecular markers in future phylogenetic considerations of Clibanarius.

4.2 Phylogenetic relationships

Asakura (2005) proposed that members of Clibanarius could be separated into two morphological groups according to the relative length of dactyls and propodi. One group contained species having stout second and third pereopods with the dactyls shorter than the propodi. The other group was composed of species with the dactyls as long as or slightly longer than the propodi, resulting in slender second and third pereopods. However, the phylogenetic results herein contradicts Asakura (2005) morphological classification scheme for Clibanarius. C. rutilus with the dactyl shorter than the propodus forms a sister relationship to C. infraspinatus with the dactyl longer than the propodus. The same pattern is found in the close relationships with C. merguiensis, C. englaucus and C. humilis, in which C. merguiensis and C. englaucus with the dactyl shorter than the propodus cluster as a sister taxon, forming a clade with C. humilis with the dactyl equal to the propodus. Overall, the present molecular data indicates that the relative length of dactyls and propodi is not a phylogenetically significant morphological character for Clibanarius.

To date, there are limited studies on the phylogeny within Clibanarius. Previous phylogenetic analyses of the Anomura and the other genera of diogenid inferred from male reproductive system or molecular markers showed low phylogenetic resolution for interspecific relationships of Clibanarius, in which Clibanarius species were regarded as a supplement (Tudge, 1997; Mantelatto et al., 2006; Tirelli et al., 2008, 2010; Tsang et al., 2011). There were only three researches exclusively explored the molecular phylogeny and morphological characteristics of Clibanarius, all which produced insufficient information on Clibanarius phylogeny and might be caused by limited taxon sampling and inadequate molecular data (Hirose et al., 2010; Negri et al., 2014; Yoshikawa et al., 2018). Our phylogenetic analysis based on all four-gene concatenated dataset presents relatively high support values for more than half of the branch nodes and the major pattern of phylogenetic relationships amongst all studied taxa is consistent across the analytic methods. The present research successfully sheds more lights on the interspecific relationships of this genus. The following interspecific relationships are drawn from analyses of all four-gene concatenated dataset: (ⅰ) C. rutilus and C. infraspinatus show a close relationship, whereas C. snelliusi is suggested to have a sister relationship with them; (ⅱ) C. longitarsus and C. striolatus are shown to be closely related and in accordance with previous studies (Hirose et al., 2010; Yoshikawa et al., 2018). Hirose et al. (2010) and Sha et al. (2015) also predicted this close affinity based on the morphological characters with the dactyl longer than the propodus as well as their habitats; (ⅲ) C. virescens clusters with the group including C. englaucus, C. merguiensis and C. humilis with high support values. The close association of the latter three species is consistent with previous molecular analyses and corroborated by their similarity in general morphology (Osawa and Chan, 2009; Hirose et al., 2010). For example, C. englaucus and C. merguiensis have common characters in dactyl, rostrum, cornea, and ocular acicle, and harbor similar inhabit environments (e.g., rocky shore, coral reef or sea grass bed Sha et al., 2015). The close affinity between C. englaucus and C. humilis was also shown in the tree topologies of Yoshikawa et al. (2018). Furthermore, we find that C. rutilus, C. infraspinatus and C. snelliusi have the telson with barely detectable median cleft, C. longitarsus and C. striolatus have the telson with very slender median cleft, while C. virescens, C. englaucus, C. merguiensis, and C. humilis have the telson with shallow median cleft (McLaughlin et al., 2007; Sha et al., 2015). Thus, we speculate that the morphological characteristic of the median cleft of the telson may be phylogenetically important for Clibanarius. In addition, the remaining two species C. eurysternus and C. corallinus have prominent U-shaped median cleft of the telson, which differs from the other Clibanarius species with shallow or/and slender median cleft. However, no agreement in the branching position of C. eurysternus and C. corallinus is obtained in the current research. The inclusion of more morphological data and extensive taxon collection are necessary to address this question in the future.

5 CONCLUSION

In this study, we investigate the phylogenetic relationships of Clibanarius based on mitochondrial and nuclear sequence data using a multigene phylogenetic approach. In contrast to Asakura's morphological classification scheme for Clibanarius, the present molecular data provides evidence that the relative length of dactyls and propodi is not phylogenetically significant for Clibanarius. Our study is the first multigene analysis on the phylogeny of Clibanarius and the phylogenetic relationships among some taxa are well resolved. We suggest that (ⅰ) C. snelliusi has a sister relationship with C. rutilus and C. infraspinatus, and the latter two taxa show closer relationships; (ⅱ) C. longitarsus and C. striolatus are shown to be closely related; (ⅲ) C. virescens is divergent from C. merguiensis, C. englaucus, and C. humilis, in which C. merguiensis and C. englaucus are sister groups with C. humilis in our research, reflecting a close relationship among them. We also speculate that the characteristic of the median cleft of the telson may be phylogenetically important for Clibanarius. In addition, our results demonstrate the importance of a large number of nucleotide sites as well as the inclusion of multiple genes with distinct evolutionary histories in the phylogenetic reconstruction of relationships among Clibanarius species. Further studies including more data from genome, extensive sampling and comparison with morphological information are necessary to elucidate phylogenetic relationships within Clibanarius.

6 DATA AVAILABILTY STATEMENT

The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request. Sequence data that support the findings of this study have been deposited in GenBank under the accession numbers MH508045–MH508087 and MK076123–MK076149.