Colour variation of the intertidal hermit crab Clibanarius virescens considering growth stage, geographic area in the Indo–West Pacific Ocean, and molecular phylogeny

Abstract Members of Clibanarius virescens show considerable intraspecific colour variation, including colouration of the second/third pereopods (green/white) and the dactyls of the second/third pereopods (with or without dark bands/patches). However, factors inducing these colour variations have not yet been elucidated. Here, we investigated the occurrence of colour variation in this species with particular emphasis on change of colouration associated with growth stage and region in specimens from tropical/subtropical to warm temperate areas in the Indo–West Pacific, including evidence from molecular phylogeny based on the cytochrome c oxidase subunit I (COI). We have, then, clarified that the colouration on the pereopod dactyls gradually changed from solid colour (yellow/white) to having dark-coloured area(s) or transverse band(s) as a result of the growth stage. The frequency of occurrence of the solid colour dactyls was higher than those of other colour types in tropical regions. Our results also indicated that the white ambulatory leg type was the colouration type that was frequently seen in juvenile stages. However, significant genetic differences were not detected between each colouration determined by molecular analysis of samples from 14 localities in the Indo-West Pacific region; in contrast, two genetically differentiated regional populations (North Australia; Phuket, Thailand; and Lombok, Indonesia) were detected. The present study, therefore, emphasizes the necessity for further study on the colour variation of marine animals focusing on growth stages and regional differences, with molecular data to facilitate the research on adaptation and/or speciation, especially in geographically widely distributed species.


Introduction
The relationship between colouration and phylogeny/natural history in organisms is an attractive topic in evolutionary biology, especially considering natural selection, adaptation and systematics (Olendorf et al., 2006;Stevens & Merilaita, 2009;Bowen et al., 2013;Green et al., 2019). In animals, body colouration has various functions, including camouflage (Stevens & Merilaita, 2009), social communication (Detto et al., 2006(Detto et al., , 2008, mimicry (Randall, 2005), and thermoregulation (Silbiger & Munguia, 2008;Kronstadt et al., 2013). It can also be an indication of reproductive isolation by serving as a cue in mate recognition (Seehausen & van Alphen, 1998). As a taxonomic tool, colouration can often be much more reliable than morphological data, when it comes to identifying the species boundaries of sister species and confirming the earliest stages of speciation in both vertebrates and invertebrates (Williams, 2000;Mathews & Anker, 2008;DiBattista et al., 2012).
Intraspecific colour differences in animal species have often been recorded, such as ontogenetic colour changes, graphical/physiological differences, and plasticity in colouration, in a broad range of animal taxa (Strong, 1975;Hill, 1999;Hanlon, 2007;Takeshita, 2019). Moreover, there are various selective pressures inducing colour morphs. In the context of pattern matching, body colouration sometimes varies considerably depending on the substrate colour, the condition of environments, and season, and function as camouflage from visual predators (Russell & Dierssen, 2018;Green et al., 2019). In terms of mating success, conspicuous colouration increases the detectability and attractivity for the counter sex, because such colouration often is originally reflected by the nutritional and/or maturity condition of individuals (Kodric-Brown, 1989;Guilford & Dawkins, 1991;Hill, 1999). Colour variation also helps in temperature regulation in various environmental conditions (Pearson, 1977;Bittner et al., 2002;Munguia et al., 2013). Intraspecific variation has led to confusion due to the lack of knowledge of the cause of the variation and the boundaries between intraspecific and interspecific differences in their systematic study (Williams, 2000;Mathews & Anker, 2008;Yoshikawa et al., 2018). Altogether, a study on intraspecific colour variation is an attractive theme in biology and helps in the understanding of the evolution, natural selection, speciation and environmental adaptation of animals. However, little is known about the intraspecific colour differences in intertidal hermit crabs.
Intraspecific variations in colouration have been documented in species of some intertidal hermit crabs. For example, the diogenid hermit crab (Calcinus morgani Rahayu & Forest, 1999) typically possesses dark brown chelipeds and ambulatory pereopods, although a specimen from a middle region of Japan exhibited purple or reddish-purple colouration on those appendages (Asakura, 2002). Moreover, five colour patterns have been identified on the second and third pereopods in Clibanarius specimens that can be assigned to the species Clibanarius merguiensis de Man, 1888 (Rahayu & Forest, 1993;Yoshikawa, personal observation). Colour difference is sometimes reflected by incipient speciation. For instance, an intertidal pagurid hermit crab (Pagurus minutus Hess, 1865) is a species that exhibits several colour variations on the lateral surfaces of its ambulatory legs which reflect genetically differentiated local populations; P. minutus are possibly separated into three groups at the species level (Jung et al., 2018).
Distinguishing between intraspecific and interspecific colour variation in hermit crabs is often difficult and has led to misidentification. For instance, Calcinus areolatus Rahayu & Forest, 1999 was later confirmed to be an early stage of C. morgani (Komai, 2004); and Clibanarius symmetricus (Randall, 1840), resurrected through molecular analysis, was previously considered to be a colour variant of C. vittatus (Bosc, 1802) (Negri et al., 2014). Following re-examination of the holotype of Clibanarius sachalinicus Kobjakova, 1955, it was revealed that C. sachalinicus was a junior synonym of Clibanarius virescens (Krauss, 1843) (Marin, 2016). Regarding the colour variants of another intertidal hermit crab (Calcinus morgani), Asakura (2002) suggested that environmental differences in its habitat could cause the colour variations in this species, as the colour-variant specimen was collected from the northern limit of this species' distribution.
Geographically widely distributed marine species often have genetically differentiated local populations or genetic structures, and the geographic differentiation pattern varies between animal taxa and even between closely related species (Paul et al., 2006;Crandall et al., 2008;Bowen et al., 2013Bowen et al., , 2016. Moreover, if the climatic conditions could be somewhat reflected in colouration, then detecting the cause of colour variation and frequency differences, especially in a species with a wide geographic range, might lead to good indicators for environmental research into global climatic changes. Knowledge of the genetic populations with the description of colour variants is, therefore, essential for conservation of local marine animal populations, as well as for understanding population differentiation and speciation. Intraspecific variation in colouration has been documented in Clibanarius virescens; this species is geographically widely distributed in tropical and warm temperate areas in the Indo-West Pacific region (Krauss, 1843;Fize & Seréne, 1955;Morgan, 1987Morgan, , 1989Haig & Ball, 1988;Forest, 1993;McLaughlin, 2002;Hirose et al., 2010;Poupin et al., 2013;Yoshikawa et al., 2018Yoshikawa et al., , 2019, and, therefore, can be suitable for discussion of geographic differences. Yoshikawa et al. (2018) reported individuals with whitish second and third pereopods (W-PCM; white propodus, carpus and merus) in this species from Nansei Island, Japan; the standard colour on the same surface is usually dark green or dark olive green. Yoshikawa et al. (2018) did not classify the W-PCM colour type as juvenile colouration because the size of the W-PCM specimens was very similar to the size of ovigerous females of normal-coloured individuals. Moreover, Morgan (1988) also reported three colouration patterns on the second and third pereopod dactyls of this species: Type A, a brown or bluish-black annular band present in the middle portion on a whitish-yellow background, which is the most common type for this species; Type B, dark dorsal and sometimes dark ventral patches present on a whitish-yellow background; and Type C, the absence of dark rings or dark patches on a whitish-yellow background (Figure 1). When we use the terminology 'dactyl' in the following section, it indicates the 'dactyls of the second and third pereopods', since we focus on the colouration of those segments in this study ( Figure 2). However, factors inducing this colour variation have not yet been investigated considering genetic structure, sex, body size, growth stages and geographic distribution.
Few descriptive studies, however, have been conducted into the factors inducing colour variations of hermit crabs, including genetic characters, sexual dimorphism, differences among the growth stages and inter-regional differences. Thus, the questions we consider in this study, using the widely geographically distributed species Clibanarius virescens, are: (1) how do local populations with different colouration in geographically widely distributed
hermit crabs differ genetically? And (2) how do the colourations differ between populations in different geographical regions?

Sample collection and examination
Field sampling of C. virescens was conducted on rocky intertidal shores of 14 Indo-West Pacific localities from April 2015 to August 2019 ( Figure 3). Individuals of C. virescens were randomly collected by hand from rocky shores one hour before/after the lowest tidal time on each sampling day. After cracking their shells, we identified the species and measured the shield length (SL) and the dactyl length (DL) of the second and third pereopods of specimens using callipers or ImageJ software 1.38 (Abramoff et al., 2004) ( Figure 2). The sex was determined by the position of the gonopores. In order to compare differences in the mature sizes of specimens of each colouration and sampling station, the shield size (SL) of ovigerous females was used as an indicator of mature body size. The colouration of the pereopods was recorded following the abovementioned three-colour types, as defined by Morgan (1988). Further, the combinations of the colour types in the dactyls of the second and third pereopods were categorized into seven groups (AA, AB, ABC, AC, BB, BC and CC; hereinafter referred to as 'dactyl-colouration-groups, DCGs'). The AA, BB, and CC represent the individuals that were composed of only A-dactyl, B-dactyl and C-dactyl type on the pereopods, respectively. The AB, ABC, AC and BC represent the individuals with each respective dactyl type.
The general colour of the pereopods was also recorded. This species generally has dark green colouration on the propodi, carpi and meri of the second/third pereopods (G-PCM), so whitish ambulatory pereopods (W-PCM) were reported as the colour variation of this species (Yoshikawa et al., 2018). Moreover, during this survey, the blue-coloured propodi, carpi and meri (B-PCM), which were previously recognized by Asakura A. (unpublished), were also recorded as a colour variation of this species (Figure 4). In summary, we investigated the size difference and occurrence between G-PCM, B-PCM and W-PCM in the present study.
In total, the ambulatory pereopods of 1209 C. virescens individuals were examined. These included 1139 and 1045 of the right second and third pereopod, and 1118 and 813 of the left second and third pereopod, respectively. Of these, 1164 individuals were used to investigate the body size difference between W-PCM and B/G-PCM pereopods colouration; this total included 237 W-PCM individuals, 99 B-PCM individuals and 828 G-PCM individuals. The pereopods colouration of 45 individuals could not be identified because these specimens had already been fixed in 99% ethanol in the field; three-colour types of dactyls could be identified even after being fixed in ethanol.
The specimens used in this study were deposited in the Seto Marine Biological Laboratory, Kyoto University, Shirahama, Japan; the Crustacean Collection of the Lee Kong Chian Natural History Museum, Singapore; the Muséum National D'histoire Naturelle, Paris, France, and the Museum of Tropical Queensland, Queensland Museum (Supplementary Table S1). Later, the newly obtained sequences were registered in the DNA data bank of Japan (DDBJ); their accession numbers are shown in Table S1 with the specimens' locality data.

Phylogenetic analysis
Sequences of the COI genes were aligned using the program Crustal W (Thompson et al., 1994) implemented in MEGA version 7 (Kumar et al., 2016) with the default settings. Poorly aligned regions were excluded from the final alignments. Phylogenetic trees were constructed using maximum likelihood (ML) methods. ML analyses were performed using RAxML (Stamatakis, 2006) as implemented in raxmlGUI 1.31 (Silvestro & Michalak, 2012). The robustness of the ML tree was evaluated with 10,000 bootstrap replications. The datasets were partitioned by codon, and the GTRGAMMA model was implemented. The tree included an outgroup using the sequence of the most closely related species, C. englaucus, following previous studies (Tsang et al., 2011;Yoshikawa et al., 2018Yoshikawa et al., , 2019.

DNA barcoding of mitochondrial COI region
Sequence data were aligned using Clustal W (Thompson et al., 1994) implemented in MEGA version 7 (Kumar et al., 2016) with the default alignment parameters and manually adjusted to minimize mismatches. The number of polymorphic sites (K), number of haplotypes (Nh), number of private haplotypes, haplotype diversity (h; Nei, 1987), and nucleotide diversity (π; Nei & Li, 1979) for each geographic area were calculated using Arlequin version 3.5 (Excoffier & Lischer, 2010).

Statistical analysis of growth size and dactyl size
For statistical analysis, a pairwise t-test with Bonferroni adjustments was conducted to visualize the differences in hermit crab body size and dactyl size between each colour type and geographic region. Additionally, Welch's t-test was conducted to elucidate the

Genetic identification and molecular analyses
Two hundred and fourteen specimens including all of the colour variations were genetically identified as C. virescens with GenBank similarities above 95% according to the resulting BLAST search. The sequences were deposited in the DDBJ with their accession numbers (Supplementary Table S1). Molecular analyses indicated that all of the DCGs sorted on the dactyl colour type (AA, AB, ABC, AC, B, BC and C) and all of the colour types of the pereopods (i.e. G-PCM, B-PCM) were included in the same clade as C. virescens ( Figure 5). A total of 164 haplotypes were found in the C. virescens populations used in this study (Supplementary  Table S1). In total, 179 variable sites and 203 mutations were found. Haplotype diversity (h) and nucleotide diversity (π) of the sampling stations ranged from 0.89 to 1.0 (mean = 0.98) and 0.01 to 0.02 (mean = 0.01), respectively. Additionally, haplotype diversity and nucleotide diversity of each phylogenetic clade ranged from 0.98 to 1.0 (mean = 0.99) and 0.01 to 0.03 (mean = 0.02), respectively. All of the molecular diversity parameters for C. virescens populations and each clade from the ML tree based on the COI region are shown in Table 2.
The median-joining network of C. virescens presented a complex haplotype network ( Figure 6) but suggested the existence of three differentiated populations in C. virescens. This result coincides with the results of ML analysis, and we categorized those as clades A, B and C, respectively (BS = 100) ( Figure 5). The central haplotype, PG-1, contained many unique haplotypes, with 18 of the 164 haplotypes sharing more than one haplotype or missing haplotype. The most common haplotype (H5) of PG-1 consisted of 16 individuals from all geographic regions (Supplementary Table S1). Closely related haplotypes generally consisted of individuals from distinct areas and did not reveal a clear geographic pattern. In contrast, 23 specimens from Phuket, Thailand, and one sample from Queensland, Australia, shared a unique haplotype distinct from the most common haplotype (PG-2, BS = 85). Moreover, three individuals from Lombok Island, Indonesia (PG-3, BS = 94), also had distinctive haplotypes from the central haplotype. Not all samples from North Australia, Indonesia and Thailand are included in PG-1, but some of them were included in PG-1. Regarding the ML tree, the haplotypes included in PG-1 were not located in the comparably close position of PG-2. Therefore, all of the haplotypes of PG-1 were essentially found in all of the sampling regions of this study. However, some unique haplotypes were found only from North Australia, Indonesia and Thailand.

Description of the haplotype relationships
Three haplotype groups were detected, and these coincided with the result of the ML phylogenetic tree. Therefore, the significant components of haplotypes in these groups are described below. Haplotype number and accession number (DDBJ) are shown in Supplementary Table S1, and their relationships were defined using the median-joining haplotype network ( Figure 6).
In PG-3, three haplotypes (H69, H72 and H73) were found, and each of these had at least one missing haplotype and two mutations. The most closely related haplotypes from PG-1 were H90 (Mie, SMBL-V0337) and H112 (Pangandaran, SMBL-V0244) with 11 mutations and one missing haplotype. We presented the haplotype composition of each sampling station in Figure 6 with classification using three PGs from the results of the phylogenetic analysis using the ML method. It also showed that the haplotype of PG-2 appeared in North Australia and Thailand.  Table 2. Results of the analysis of molecular diversity parameters for Clibanarius virescens populations and each clade from the ML tree based on COI region: number of sequences analysed (N), variable sits (S), total number of mutations (Eta), average number of nucleotide differences (k), number of haplotypes (h), haplotype diversity (Hd), haplotype diversity standard deviation (SD Hd), nucleotide diversity (π), nucleotide diversity standard deviation (SDπ), Fu's Fs statistic (Fu's Fs)  Tajima's D values and Fu's Fs indicated negative but nonsignificant values, which suggests the existence of many silent mutations (i.e. mutations that do not change the amino acid sequence), thus producing low heterozygosity as there were some synonymous mutations (Table 2). These negative values agreed with the high DNA sequence diversity and haplotype diversity of C. virescens ( Figure 6; Table 2; Supplementary  Table S1). Moreover, it generally indicated genetically mixed populations in which mutations were accumulated over long periods and, thus, allowed for the maintenance of these high levels of diversity. The presence of the most abundant haplotypes in PG-1, together with the high nucleotide diversity (π) in each station, suggested high gene flow and null population differentiation of the populations from the West Pacific Ocean.

The comparison of the DL among three colourations
The frequencies for each colouration are shown in Figures 7  and 8. The colourations appear in the order C-type, B-type and A-type from small to large crabs. The results of the Kruskal-Wallis χ 2 test indicated that there were significant differences in DL between colourations in both sexes (Table 3). We then performed a pairwise t-test to visualize the size differences between A-type, B-type and C-type dactyls. Significant size differences were found in almost all of the combinations of DL, except the A-type and B-type of the females' left third pereopods (Table 3). These results also showed that the frequency of the three-colour variations of the dactyls depended on the growth stage. All juveniles (SL: 1.5-2.0 mm) had the C-type pattern. After reaching specific adult sizes (∼2.5 mm in females and 2.0 mm in males), A-type and B-type patterns emerge. While the frequency of the C-type pattern decreased and the A-type increased with growth in both sexes, the rate of the B-type pattern remained at ∼20-40%; the DL was in the range of 2.0-5.5 mm in females and 2.5-7.0 mm in males.

The body size differences among the seven DCG
Examined specimens included 275 specimens of AA, 111 of AB, 6 of ABC, 13 of AC, 201 of BB, 52 of BC and 551 of CC. The results of the Kruskal-Wallis χ 2 test indicated that SL significantly differed among the seven categories AA, AB, ABC, AC, BB, BC and CC; and between females and males (Table 4). We performed a pairwise t-test to visualize the size differences among these colour categories, and the results are shown in Table 4. Significant differences in female body size were found between AA-BB, AA-BC, AA-CC, AB-BC, AB-CC and BB-BC (P ≤ 0.01). In males, significant differences in male body size were found, except between AB-AC, AC-BB, AA-BC and BB-BC.

The body size difference and frequency of the seven colour DCGs among each region
The differences in the frequency of the seven colourations are shown in Figures 9 (female) and 10 (male). The percentages of the individuals that possess the B and/or C dactyls were higher in tropical regions than temperate regions in both sexes (female, 72% in temperate regions (A, B and C in Figure 3) and 100% in tropical regions (D, E and F in Figure 3); male, 66% in temperate regions and 99% in tropical regions) (Figures 9 and 10). Also, the body size tended to be relatively small in the tropical population (female, temperate regions = 3.36 ± 0.69 mm, tropical regions = 2.7 ± 0.52 mm; male, temperate regions = 3.82 ± 1.56 mm, tropical regions = 3.52 ± 0.86 mm) (Figures 9-11). To understand the differences in the mature sizes between each sampling station, the SLs of 189 ovigerous females were examined (Table 5). The minimum size and mean size of the SLs from tropical regions were smaller than in the temperate regions (Tropical regions: North Australia, minimum size = 2.0 mm, mean ± SD = 2.75 ± 0.74 mm; Indonesia, minimum size = 2.1 mm, mean ± SD = 2.8 ± 0.36 mm; Thailand, minimum size = 2.07 mm, mean ± SD = 2.68 ± 0.24 mm. Temperate regions: Honshu, minimum size = 2.56 mm, mean ± SD = 3.84 ± 0.74 mm; Kyushu, minimum size = 3.12 mm, Fig. 6. The haplotype network of C. virescens. The colour of circles represents each ocean area (green = Japanese water; red = North Australia (Queensland); yellow = Indonesia (Jawa and Lombok Island); and navy blue = Thailand (Phuket)).   mean ± SD = 4.04 ± 0.71 mm; Nansei Island of Japan, minimum size = 3.39 mm, mean ± SD = 3.79 ± 0.31 mm) ( Figure 11; Table 5). The results of the Kruskal-Wallis χ 2 test on the SL of ovigerous females also indicated that SL significantly differed among the sampling stations (Table 5). The pairwise t-test also indicated that the shield size of ovigerous females was significantly different between A-E, A-F, B-D, B-E, B-F, C-D, C-E and C-F sampling stations (Table 5). Therefore, mature females in the tropical population were smaller than those in the temperate region.

Body size difference between W-PCM and B/G-PCM pereopods colouration
We investigated the size differences of B-PCM, G-PCM and W-PCM (Figure 12). A total of 1164 C. virescens individuals, which Numbers next to the colour type indicate the number of the dactyls of C. virescens investigated in the study. Double asterisks (**) and triple asterisks (***) indicate P ≤ 0.01 and P ≤ 0.001, respectively. included 591 females and 573 males, were used. This total also included 99 B-PCM individuals, 828 G-PCM individuals and 237 W-PCM individuals (Table 6). Significant bias was detected by the Kruskal-Wallis χ 2 test that was conducted on the SL between the three colourations of females and males, respectively. Significant differences in female body size were found between B/ G-PCM, G/W-PCM and G/W-PCM. In males, significant differences in male body size were found between colourations, except   for the comparison of B/W-PCM and G/W-PCM (Table 6). However, no significant difference was detected between male B/G-PCM patterns. Although the B-PCM and G-PCM colour types are mainly found in body sizes that are larger than 2.0 mm, the W-PCM type was only found between the range of ∼1.5-2.5 mm ( Figure 12). Our results indicated clear evidence which confirms that the W-PCM type is a juvenile colour type and, after reaching 2.5 mm, the colour changes to B-PCM or G-PCM (Figure 12).
Regional difference on the body size and frequency of B/G-PCM Regarding the regional difference in the frequency of B-PCM and G-PCM, since B-PCM type were not collected in Indonesian and North Australian regions during this study, the specimens from four areas, including Honshu, Shikoku, Nansei Islands, Japan, Indonesia and Thailand, were examined. First, we compared the size differences of each sampling site using t-tests. Significant size differences were detected between B-PCM and G-PCM of the female population from Thailand (P ≤ 0.05). However, no significance was detected in other populations of either sex.
To understand the relationships between the frequencies of the B/G-PCM and sampling regions, we used Fisher's Exact Test (B-PCM or G-PCM × Number of the sampling station). The results of Fisher's Exact Test (2 × 4) indicate that there is a significant bias in the occurrence of B/G and sampling regions in both sexes (female, P ≤ 0.001; male, P ≤ 0.001) (Figures 13 and 14). Therefore, considering our observations and the statistical analysis results, the frequency of the occurrence of B-PCM pereopods is higher in Phuket, Thailand.

Differences in the occurrence of dactyl colour variations
Our results indicated that the frequency of the colour variation changed from C-type to A-and/or B-type depending on the growth stage. However, because the dactyls also exhibited varied colourations even in one specimen (categorized as AB, ABC, AC and BC) (Figures 7 and 8), the colour transformation of dactyls might depend on not only the growth stage but also the stages of dactyl after autotomy. Our assessment of the seven types (AA, AB, ABC, AC, BB, BC and CC) based on their dactyl colouration indicated that the occurrence of the AA-type was different depending on the geographic area. We, therefore, concluded that there were regional differences in the frequency of the colour variations of the dactyls. This regional difference might be related to the difference in the growth/matured rate between tropical and template regions since the mean maturity size of the population inhabiting the tropical areas was significantly smaller than that of the temperate region ( Figure 11 and Table 5).
These growth/maturity patterns might follow the temperaturesize rule, which has been empirically supported in many ectotherms (Kingsolver & Huey, 2008;Forster et al., 2012). In crustacean species, the growth rate, mean size and maturity size are also different depending on the optimal environmental temperature of the species. For example, the mean size and maturity size of the fiddler crab Tubuca arcuata (De Haan, 1835) is smaller in tropical/subtropical areas than in temperate areas (Aoki et al., 2010). As another example, the sandy beach isopod Excirolana braziliensis Richardson, 1912 experiences a decrease in lifespan and an increase in natural mortality from temperate to subtropical beaches, but increases their mass at size (length-mass relationship) from subtropical to temperate beaches (Cardoso & Defeo, 2004). The study hypothesized that isopods in temperate regions (at higher latitudes) achieve larger sizes by delaying maturity; therefore, they may be able to invest more resources in reproduction and body growth than individuals in tropical water. In   C. virescens, the growth and maturity rate of individuals inhabiting tropical regions may be faster than those inhabiting temperate regions (as shown in Figure 11), and populations in tropical regions might have a higher mortality rate or a shorter lifespan. Moreover, the frequency of the A-and B-dactyl types, which were mostly exhibited on larger individuals, were lower in tropical regions. Therefore, we suggest that the transformation of dactyl colour depends entirely on their body size and not on their maturity. It is known that under global climatic changes, geographic range of species are gradually moving towards higher latitudes. The land hermit crab Coenobita purpureus Stimpson, 1858, has been shown to increase its northern limits on the Pacific coast of Japan (Sanda et al., 2019). It is possible that the frequency of the colour variation of C. virescens may follow the same pattern. Under the effect of global warming, the colour composition in southern locations might be gradually detected in more northern areas. The geographic difference in the colour morph could be the indicator for global climatic changes and, therefore, further continuous time-series monitoring will be needed. Another possibility for this regional difference of dactyl colouration is the restricted growth size of C. virescens through limited shell resources, due to the high-pressure interspecific competition for snail shells. The potential preference of shell use pattern (species-specific shell preferences) varies between hermit crab species (Reese, 1969;Imazu & Asakura, 1996), and the shell usage patterns tend to restrict growth rate and reproductive potential. Interspecific competition seemed relatively low in temperate regions, since, in Clibanarius species, only one species (C. virescens) has been reported in rocky shores in the temperate regions of Japan (Honshu, Shikoku and Kyushu) (Imazu & Asakura, 1996 and in this study). In contrast, more than two species of Clibanarius sympatrically inhabit the tropical area, in general. Therefore, interspecific competition may be higher in tropical regions. In this study, six Clibanarius species (C. corallinus, C. englaucus, C. eurysternus, C. humilis, C. snelliusi and C. striolatus) were collected in tropical regions of the northern hemisphere Pacific ocean (e.g. Amami-Oshima and Okinawa Island). One species, C. taeniatus was collected at the sampling station of Townsville, Queensland, Australia. Six species, C. corallinus, C. englaucus, C. eurysternus, C. humilis, C. merguiensis and C. striolatus, were collected from the sampling stations of Jawa and Lombok Island of Indonesia. And two species, C. cruentatus and C. merguiensis, were collected from Laem Panwa, Phuket, Thailand. Since competitive pressure might be higher in tropical regions, the growth ability of C. virescens may be restricted. Thus, the regional differences in the frequency of dactyl colourations may have occurred due to limited shell resources as a result of high interspecific competition for snail shells.
Differences in the occurrence of colour variations of pereopods Yoshikawa et al. (2018) suggested that W-PCM was not a juvenile colour. However, no quantitative investigation was conducted in Yoshikawa et al. (2018). Thus, it was not clear whether the W-PCM type was a juvenile colour type or whether the colour type of C. virescens changes depending on its growth stage. In this study, our results confirmed that this colour type is frequently seen in the juvenile stage, as it was only found with an SL size range of 1.0-2.5 mm. The colouration changes to green or blue from white as the individual grows. As Yoshikawa et al. (2018) mentioned, ovigerous females were also found close to this size (2.0-2.5 mm) (Figure 11), and the typical colouration (dark olive green) was also found for individuals with an SL size of 1.0-2.5 mm; no W-PCM individuals were found with more than a 3.0 mm SL size for both sexes ( Figure 12). Also, ovigerous females with W-PCM colouration were not found in this study, and all ovigerous females exhibited either B-PCM or G-PCM colouration. Therefore, we conclude that the W-PCM type is a juvenile C. virescens colour, and the typical colours for this species are usually found in those individuals with ∼2.5 mm in SL size.
Body colour is sometimes different between juveniles and adults of other hermit crab species, for example, Rahayu & Forest (1999) described a new species, Ca. areolatus, which superficially resembles Ca. morgani Rahayu & Forest, 1999, but with small dots on the dactyls and propodi of the second and third pereopods. However, Komai (2004) revealed that Ca. areolatus is a juvenile stage of Ca. morgani. Because the distinctively different juvenile colourations can cause hermit crab taxonomic confusion, quantitative investigation of the colour transition of hermit crabs is very important for future taxonomic and ecological studies.
In this study, we also observed a B-PCM colouration on the pereopods, which was previously recognized by Asakura A. (unpublished) (Figure 4). Clibanarius virescens was described based on specimens obtained from the coast of Natal, South Africa (Krauss, 1843), and the colouration was originally described by Krauss (1843) as 'dark olive-green'. Crustacean body colour sometimes differs in the stage immediately following
moulting (Kurup, 1964;Reaka, 1975). In some Anomura species, the carapace and chelipeds are completely hardened after ∼20 and 7 days, respectively, following moulting (Kurup, 1964). Moreover, it is also known that colouration develops best during the postmoulting stages in other crustaceans, such as stomatopod species (Reaka, 1975). However, in C. virescens, it seemed as though the colour difference does not only depend on the moulting stages. If the colour type of C. virescens did indeed depend on those stages, the B-PCM ought to have been found at every SL size at a constant rate. Although the occurrence of the B-PCM type was very low for both sexes in Japanese and Indonesian waters, the frequency of this pattern was relatively high in Phuket, Thailand (Figures 13  and 14). Excluding females, no significant SL size difference was found between green-pereopod and B-PCM individuals from Phuket, Thailand; only females with B-PCM were significantly smaller than the green-pereopod individuals ( Figure 13B; Table 6). As a result, it was difficult to conclude whether the B-PCM colouration was caused by growth stages or not. Therefore, we could not ascertain the reason for this geological difference or the significance of this colour variation. This colour type might be related to the feeding habits of the Phuket, Thailand population. Clibanarius virescens is known to climb out of tidal pools onto rock surfaces exposed to air during low tide; it experienced a peak at ∼3.0-4.0 mm shield size (Yoshikawa et al., 2020). Thus, it is also possible that the colour formation might change depending on the intensity of sunlight in each environment. Otherwise, negative frequency-dependent selection may have occurred on B-PCM and G-PCM to maintain the genetic variation in their population. Therefore, the frequency of these two colourations may be dynamically changing on a long timescale.

Genetic differences linked to colouration
Our molecular analyses suggest that the monophyly of all colour DCGs sorted following the dactyl's colour type (AA, AB, ABC, AC, B, BC and CC) and colour types G-PCM and B-PCM were included in the same clade (PG-1) of C. virescens. Although PG-2 and 3 are composed of only CC specimens, it may be biased since almost all specimens used in molecular analysis exhibited BB, BC or CC colour types. Therefore, we concluded that all of these DCGs and the colour variations of the pereopods (G-PCM and B-PCM) are the result of intraspecific variations of C. virescens, but not interspecific differences. Although no genetic differences between the G-PCM, B-PCM and each DCG were detected, the ML analyses and haplotypes network analyses showed three genetically differentiated populations (PG-1, 2 and 3) (Figures 5 and 6). PG-2 was composed of some North Australian and Thailand samples, and PG-3 was only composed of some of the samples from Lombok Island, Indonesia.

DNA barcoding and phylogeographic pattern
Clibanarius virescens was originally described based on specimens obtained from the coast of Natal, South Africa (Krauss, 1843), and has been documented on many rocky shores and coral reef localities of the Indo-West Pacific region, such as Mayotte Island (Poupin et al., 2013), Indonesia (Rahayu & Forest, 1993), Vietnam (Fize & Seréne, 1955), Thailand (McLaughlin, 2002), Guam (Haig & Ball, 1988), Somalia (Lewinsohn, 1982), Australia (Morgan, 1987(Morgan, , 1989 and Japan (Hirose et al., 2010;Yoshikawa et al., 2018Yoshikawa et al., , 2019. Although C. virescens has been recognized as one of the most geographically widely distributed hermit crabs, its geographic differentiation pattern has not yet been elucidated using molecular methods. Further, the demographic history and existence of cryptic species or populations have not been discussed yet. In our molecular analysis, PG-1 consisted of the samples from all localities, and these haplotypes were connected by some mutations with or without several missing haplotypes and showed a more complex haplotype network than other PGs. As a factor of speciation and population deviation, the sea-level fluctuation during the Pleistocene era that made the land between the Indian and the Pacific Ocean is a well-known factor for many organisms of tropical marine biota (Smith & Sandwell, 1997;Williams & Benzie, 1998;Voris, 2000). Therefore, C. virescens might have also experienced the bottleneck effect of removing the ancestral polymorphisms of the COI gene around Phuket, Thailand.
Moreover, according to our ML phylogenetic reconstruction, the haplotype from Thailand, which was in PG-1, was not located in the nearest position to PG-2 (purple branches of Figure 5). This result indicated that the Thailand haplotypes in PG-1 are not closely related to PG-2. Since some haplotypes that belonged to PG-1 (major haplotype) were also detected in Phuket, Thailand, a second genetic inflow to Phuket from the West Pacific regions might have occurred during the planktonic larvae stage by using ocean currents or by being introduced via a ship's ballast water. One North Australian haplotype that belonged to PG-2 might support this hypothesis. Endemic haplotypes (PG-3) were also found in Lombok Island, Indonesia. Near Lombok Island, the Malacca Strait is a geographic barrier that separates the genetic population of tropical marine biota on Jawa-Bali Island and Lombok Island (Williams & Benzie, 1998;Reid et al., 2006). Therefore, it may also work as a geographic barrier to induce the deviation of C. virescens populations. Otherwise, if PG-2 and PG-3 represent the early stages of speciation in Phuket, Thailand and Lombok Island, it would be difficult to see how they are reproductively isolated from the PG-1 individuals, although they occur sympatrically. Behavioural experiments on mate choice and/or more detailed population genetic analysis in these two fields, respectively, may answer the above question.
Another possibility is that the actual population might have more haplotypes connecting with the haplotypes of PG-2 and PG-3; the haplotypes belonging to PG-2 and -3 may also have been distributed between the Coral Sea to the East China Sea. The core-periphery hypothesis (da Cunha et al., 1950;Brussard, 1984) and the abundant-centre (Antonovics, 1976;Hengeveld & Heack, 1982) hypothesis are well known as expected outcomes of high genetic diversity in a geographically widely distributed cosmopolitan species. There is still the possibility that there is more genetic diversity in tropical C. virescens populations, as our molecular population analysis shows. To confirm this possibility, both the analysis of other DNA regions (e.g. several nucleus DNA markers in addition to COI) and comprehensive sampling between each survey region is necessary.

Summary and conclusion of the present study
Our research focused on the colour variants of the intertidal hermit crab C. virescens, and we investigated the factors and genetic characters impacting on colouration. The investigation was conducted with particular emphasis on speciation (e.g. detecting cryptic species/populations), phylogenetic/genetic relationships, growth stages, sex and regional differences in specimens from tropical/subtropical to warm temperate areas in the Indo-West Pacific. We, then, confirmed that the colouration on the pereopod dactyls gradually changed from a solid colour (yellow or white) to having dark-coloured area(s) or transverse band(s) as a result of the growth stage. The frequency of occurrence of the solid colour dactyls was higher than those of other colour types in tropical regions. Moreover, it was also confirmed that the white ambulatory leg type was a part of the colouration frequently seen in the juvenile stage. None of these colour variants were genetically differentiated in our phylogenetic analysis with the mitochondrial cytochrome oxidase subunit I (COI) gene. In summary, the latitudinal/growth-stage difference in the occurrence of colour variation was found in C. virescens, without a genetic distinction; the genetic structure was detected in the geographic population, not on each colour morph. To the best of our knowledge, this study is the first investigation of the factors influencing the colour variation of C. virescens, and it emphasizes the necessity for further investigations on colour variations of morphologically similar marine species linked with taxonomic and ecological studies.
Supplementary material. The supplementary material for this article can be found at https://doi.org/10.1017/S002531542000106X