Constrained body shape among highly genetically divergent allopatric lineages of the supralittoral isopod Ligia occidentalis (Oniscidea)

Abstract Multiple highly divergent lineages have been identified within Ligia occidentalis sensu lato, a rocky supralittoral isopod distributed along a ~3000 km latitudinal gradient that encompasses several proposed marine biogeographic provinces and ecoregions in the eastern Pacific. Highly divergent lineages have nonoverlapping geographic distributions, with distributional limits that generally correspond with sharp environmental changes. Crossbreeding experiments suggest postmating reproductive barriers exist among some of them, and surveys of mitochondrial and nuclear gene markers do not show evidence of hybridization. Populations are highly isolated, some of which appear to be very small; thus, the effects of drift are expected to reduce the efficiency of selection. Large genetic divergences among lineages, marked environmental differences in their ranges, reproductive isolation, and/or high isolation of populations may have resulted in morphological differences in L. occidentalis, not detected yet by traditional taxonomy. We used landmark‐based geometric morphometric analyses to test for differences in body shape among highly divergent lineages of L. occidentalis, and among populations within these lineages. We analyzed a total of 492 individuals from 53 coastal localities from the southern California Bight to Central Mexico, including the Gulf of California. We conducted discriminant function analyses (DFAs) on body shape morphometrics to assess morphological variation among genetically differentiated lineages and their populations. We also tested for associations between phylogeny and morphological variation, and whether genetic divergence is correlated to multivariate morphological divergence. We detected significant differences in body shape among highly divergent lineages, and among populations within these lineages. Nonetheless, neither lineages nor populations can be discriminated on the basis of body shape, because correct classification rates of cross‐validated DFAs were low. Genetic distance and phylogeny had weak to no effect on body shape variation. The supralittoral environment appears to exert strong stabilizing selection and/or strong functional constraints on body shape in L. occidentalis, thereby leading to morphological stasis in this isopod.


Introduction
Morphological stasis, the lack of change in gross external anatomy over long periods, is conspicuous in the fossil record and among extant taxa (Gould and Eldredge 1977;Wake et al. 1983;Bickford et al. 2007). Considered one of the most challenging problems in evolutionary biology, this phenomenon is central to understanding the gradual-ism versus saltation discrepancy (Gould and Eldredge 1977;Charlesworth et al. 1982;Wake et al. 1983;Futuyma 2009). Cryptic diversity, which is pervasive in nature (Bickford et al. 2007;Trontelj and Fiser 2009), provides remarkable opportunities to study the extent and underlying causes of morphological stasis. Studies of extant cryptic species have revealed morphological stasis among lineages that appear to have diverged long ago (e.g., up to tens of millions of years) and occupy distant geographic areas (Lee and Frost 2002;Lavoue et al. 2011). Nonetheless, in certain taxa formerly regarded as examples of morphological stasis, the application of geometric morphometric approaches uncovered previously unknown and lineage-diagnostic divergence in external morphology (Strumbauer and Meyer 1992;Maderbacher et al. 2008;van Steenberge et al. 2015). Therefore, variation in gross external morphology must be carefully evaluated before a taxon is deemed "morphologically static." Members of the semiterrestrial isopod genus Ligia appear to have experienced morphological stasis. DNA sequence data have revealed high levels of isolation and genetic allopatric differentiation within recognized coastal species of this genus from separate regions (i.e., Hawaiian archipelago, Gulf of California-adjacent areas, Caribbean, and Mediterranean), each of which is suggested to harbor a cryptic species complex (Hurtado et al. 2010;Santamaria et al. 2013Santamaria et al. , 2014Hurtado et al. unpublished). Traits such as direct development (Carefoot and Taylor 1995), which restrict coastal species of Ligia to a vertically narrow strip of the coast (i.e., rocky upper intertidal and supralittoral), have likely contributed to these high levels of allopatric genetic differentiation (Hurtado et al. 2010). Members of Ligia are among the few animals that have adapted to live exclusively within this environment, which is characterized by extreme conditions . These include regular exposure to a broad range of temperatures, air humidity, and water salinity levels (from rain, wave splash, storm surge, and tides), and to predation by terrestrial, aerial, and marine animals (Menge 1976;Brown 2001;Ellis et al. 2007;Castilla et al. 2008;Donahue et al. 2009). Such conditions might impose strong stabilizing selection on morphology, thereby limiting the morphological divergence that can accompany genetic differentiation.
The rocky intertidal isopod Ligia occidentalis (Dana 1853) represents the most striking example of morphologically cryptic diversity reported to date within the genus (Fig. 1). It is distributed on the Pacific coast of North America from southern Oregon to Central Mexico, including the Gulf of California. This isopod is currently recognized as a single species for which no junior synonyms have been proposed (Espinosa-P erez and Hendrickx 2001;Schmalfuss 2003). Phylogeographic analyses of L. occidentalis, however, revealed the existence of highly divergent lineages that likely represent a cryptic species complex (hereafter L. occidentalis sensu lato). Divergences among lineages are as high as 29.9% (Kimura-2-parameter; K2P) for the Cytochrome Oxidase I gene (COI), where > 60% of pairwise comparisons among localities exhibit COI K2P divergences > 10% (Hurtado et al. 2010). Deep divergences among lineages suggest a long evolutionary history of L. occidentalis s. l. in the region, possibly since the Miocene, whereas high levels of genetic differentiation among populations indicate gene flow is severely restricted, even over small geographical distances, implying long-standing isolation of populations (Hurtado et al. 2010). Crossbreeding experiments suggest postmating reproductive barriers may exist among some lineages (McGill 1978), and surveys of mitochondrial and nuclear gene markers do not show evidence of hybridization among divergent lineages (Eberl et al. 2013;Hurtado et al., unpublished data).
Divergent lineages of L. occidentalis s. l. are exposed to markedly different environmental conditions, which could entail natural selection promoting morphological divergence. Ligia occidentalis s. l. is distributed along ã 3000 km latitudinal gradient that encompasses several proposed marine biogeographic provinces and ecoregions in the eastern Pacific (Spalding et al. 2007;Robertson and Cramer 2009;Briggs and Bowen 2012). In general, the eight main lineages have nonoverlapping geographic distributions, with distributional limits that generally correspond with sharp environmental changes. The geographical limit between the two most divergent Pacific clades of L. occidentalis s. l. (20-25% divergence for COI) distributed from southern Oregon to the southern Baja California Peninsula occurs at the Point Conception biogeographical boundary, which separates the Oregonian and Californian marine zoogeographical provinces (Eberl et al. 2013). The geographical limit between these two main clades largely reflects the changes in sea surface temperature that define the Point Conception biogeographical boundary along the shores of both the mainland and the Northern California Channel Islands (Eberl et al. 2013).
Within the Gulf of California, sharp environmental differences are also observed in the distribution of the two most divergent lineages of Ligia (15-26% K2P COI) occupying this basin (Hurtado et al. 2010). One lineage occurs in the Northern Gulf of California (from the mouth of the Colorado River to the Midriff islands), whereas the other is distributed in the Southern Gulf of California. The Northern Gulf is characterized by strong seasonal variation in water temperatures (>30°C in the summer; 8-12°C in the winter), large tidal regimes (up to 10 m), and high summertime salinities (35-40 ppt), whereas the Southern Gulf is characterized by somewhat lower salinities, smaller tidal regimes, and moderate seasonal variation in water temperatures (30-32°C in the summer; 18-20°C in the winter) (Brusca 2007 and references therein).
Large genetic divergences among lineages, reproductive isolation among them, high isolation of populations, and/ or marked environmental differences in their ranges may have promoted morphological differentiation in L. occidentalis s. l. not detected yet by traditional taxonomy. Use of more sophisticated approaches, such as geometric morphometrics, may detect differences among lineages of this isopod. Geometric morphometric analyses have been useful for discriminating cryptic lineages of other crustaceans (Silva et al. 2010a,b;Zuykova et al. 2012), as well as in an array of other animal taxa (e.g., Carvajal-Rodr ıguez et al. 2006;Francuski et al. 2009;Milankov et al. 2009;Mitrovski-Bogdanovic et al. 2013;Schmieder et al. 2015). Remarkable body shape differences can be attained rapidly in response to environmental variation (e.g., the freshwater isopod Asellus aquaticus; Eroukhmanoff and Svensson 2009). Body shape is one of the most important and comprehensive features of an organism's phenotype (Ingram 2015), and is relevant to ecological traits of isopods (e.g., Schmalfuss 1984;Broly et al. 2015). The presence of a rigid exoskeleton in crustaceans facilitates unambiguous placement of landmarks for geometric morphometrics analyses of body shape.
In this study, we used landmark-based geometric morphometric analyses to test for differences in body shape among highly divergent lineages of L. occidentalis s. l., and among populations within these lineages. We used landmarks that captured taxonomically informative regions that have been used to distinguish Ligia species. We conducted discriminant function analyses to test whether body shape morphometrics can be used to diagnose genetically differentiated lineages of L. occidentalis s. l. We used thin-plate-spline transformations to describe the general shapes of individuals within lineages and make comparisons among lineages. We also tested (1) for associations between phylogeny and morphological variation, and (2) whether genetic divergence is correlated to multivariate morphological divergence. Our study contributes to understanding the constraints on morphological evolution in a cosmopolitan genus characterized by high levels of cryptic diversity and an extreme habitat (i.e., the supralittoral).

Samples
We analyzed a total of 492 Ligia individuals from 53 Pacific localities distributed between central California and central Mexico, including the Gulf of California ( Fig. 2; Table 1). Samples were collected by hand during 2007-2010 and stored in 100% ethanol under permits from California Department of Fish and Game (USA) No. 9881, and Comisi on Nacional de Acuacultura y Pesca (Mexico) No. DGOPA.l0337.020908.2952. Once in the laboratory, they were stored in ethanol at À80°C until dissection. These specimens were part of the original collections obtained for the Ligia phylogeographic study of this region by Hurtado et al. (2010) and represent the eight main highly divergent lineages (clades) identified in that study, which, for the most part, have nonoverlapping distributions (Hurtado et al. 2010;Eberl et al. 2013). For consistency among studies, we use the same names for these lineages as in Hurtado et al. (2010). (1)
We characterized body shape by digitizing 27 landmarks, using TpsDig v2.16 (Rohlf 2004), on the periphery of Ligia bodies (Fig. 3). Care was taken to include landmarks that captured taxonomically informative regions and that can be measured unambiguously. Landmarks were placed on medial and lateral boundaries of the eyes at the body periphery. These landmarks capture the relative size of the eyes and the distance between them, both characters used to distinguish Ligia species (Taiti et al. 2003). Landmarks were also placed on lateral posterior tergite tips to aid in characterizing relative width of each body segment and overall body shape, also important in Ligia taxonomy (Jackson 1922;Schultz and Johnson 1984;Lee 1994;Taiti et al. 2003; Khalaji-Pirbalouty and W€ agele 2010). Finally, landmarks were placed at the posterior tip and the lateral posterior points of the pleotelson. Relationships between these landmarks capture the shape of the pleotelson, another trait used in Ligia taxonomy (Schultz 1974;Taiti et al. 2003; Khalaji-Pirbalouty and W€ agele 2010).
As the body plan of Ligia is bilaterally symmetric, all but the pleotelson tip landmarks are anatomically homologous and should not be treated as independent in statistical analyses. We therefore reflected and averaged homologous landmarks across the midline (Zelditch et al. 2004), which was defined as a line connecting the pleotelson tip and the midpoint between the medial eye landmarks. Corrected landmarks were centered, scaled, and rotated, to best align with the consensus (i.e., average landmark configuration), using the method of generalized least squares, and projected to a flat (i.e., tangent) shape space using tpsRelw v1.49 (Rohlf 2006). We calculated principal components (PCs) of aligned coordinates to yield orthogonal shape variables. We also estimated logcorrected centroid sizes (summed squared distances of landmarks from the centroid; Bookstein 1991) to use as a measure of body size.

Statistical analyses
We conducted full factorial MANCOVA analyses of shape variables as a function of Lineage (i.e., Clade), Population nested within Lineage, Sex, Size, and all interactions, to discern the meaningful correlates of body shape. When interaction terms were not significant, we removed them from the model, hierarchically by order (i.e., from more complex to simpler), and repeated the analyses (Engqvist 2005). We estimated effect strengths for all terms in the final model by calculating partial eta-squared values (g 2 P ), which is the multivariate analog of R² in simple regression models (Tabachnick et al. 2001).
We explored differences among lineages, and among populations within lineages with discriminant function analyses (DFAs). To focus exclusively on between group differences, we first accounted for other predictors by conducting a preliminary MANCOVA based on Size and Sex and saving residual variation. Residuals were used in DFAs to focus discrimination solely on Lineage effects (Langerhans and DeWitt 2004). We tested whether all groups in our data shared a covariance matrix with the Box's M test. Although quadratic DFAs do not assume a homogeneous covariance matrix, singularities in the data matrix may prevent their correct use. In cases where neither linear nor quadratic DFAs could be correctly applied, we used regularized DFAs (Friedman 1989), as a compromise approach. We determined the best combination of k (i.e., the degree of shrinkage of the individual class covariance matrix estimates toward the pooled estimate) and c (i.e., the degree of shrinkage toward a multiple of the identity matrix) values by evaluating the risk of misclassification under several combinations of these parameters as suggested by Friedman (1989). We specified equal prior probabilities for each lineage in all DFAs analyses. Using this procedure, we attempted to assign individuals to (1) their clade of origin and (2) their population of origin within their corresponding clade. All results were validated using leave-one-out cross-validation (LOOCV).
We tested for associations between phylogeny and morphological variation by estimating Pagel's k (Pagel 1999) and Blomberg's K (Blomberg et al. 2003) for all shape variables, using Ives et al.'s (2007) method to account for multiple observations per terminal branch. The use of shape PCs in these analyses is justified, as they are aligned to the main axis of variation and maintain interobject distances (Perez et al. 2011). Both statistics provide a univariate measure of the strength of phylogenetic signal in the data, with values close to zero indicating no phylogenetic signal, and values close to one indicating the character has evolved under a Brownian motion (BM, i.e., phylogenetic signal explains the observed patterns). We tested whether observed k values were statistically different from those expected under a null model (i.e., BM = 0) and a fully Brownian model (i.e., BM = 1) using a likelihood ratio test. In addition, we tested whether observed K values departed from the null hypothesis of no phylogenetic signal using 10,000 permutations (Blomberg et al. 2003). All computations were carried out using the Picante (Kembel et al. 2010) and GEIGER (Harmon et al. 2008) packages in R.
We also tested whether genetic divergence is related to multivariate morphological divergence. We estimated genetic distances from COI sequences published by Hurtado et al. (2010) using the K2P model in PAUP* (Swofford 2003). We calculated pairwise Euclidean distances for all localities on the residual variation (see above) using PopTools v. 3.2 (Hood 2010) in Microsoft â Excel. We tested for correlations between COI K2P distances and Euclidean Morphological Distances (i.e., a Mantel test) and estimated P-values by permutation. All statistical tests were carried out in JMP v9.0.1 (SAS Institute Inc., Cary, NC).
To visualize shape differences among lineages, we produced thin-plate-spline transformations (average shape deformation for each lineage from the consensus shape) of LM positions in tpsRegr v1.37 (Rohlf 2005) by entering our MANCOVA design matrix as the independent variable and our symmetrical landmark constellations as the dependent variables. We also produced transformations for the smallest and largest individuals for each clade and in the dataset overall to visualize shape differences between specimens of differing body size. We used these transformations to describe the general shapes of individuals within lineages and make comparisons among lineages.

MANCOVA
Principal components analysis generated 24 nonzero eigenvectors. The first eleven PCs accounted for 95.4% of the variance and were included in subsequent analyses, whereas the other thirteen were discarded. The full factorial MANCOVA yielded no significant three-or four-way interaction terms, which were removed prior to repeating the analysis. This simpler MANCOVA model (  (Table 2). Of these, the only effects with a partial eta square (g 2 P ) value above 0.2 were Lineage, Size, and the interaction term Population 9 Size [Lineage] ( Table 2). The main effect of Sex was not significant (Table 2).

Discriminant function analyses
Results of the Box's M test indicated that covariance matrices are heterogeneous across lineages (Box's M: 1123.2, df error = 72611.9, P < 0.0001), suggesting linear DFAs were inappropriate for our dataset. We could not implement quadratic DFAs, however, due to singularities in our data matrix. Therefore, we implemented regularized DFAs. We used k and c values of 0.1 in the final analyses, as low values for these parameters are recommended when covariances are different, data are abundant, and when variables may be correlated (SAS Institute Inc. 2010). Also, these values produced the lowest misclassification rates in preliminary analyses under a variety of k and c combinations, another criterion for parameter selection (Friedman 1989). Regularized DFAs of residuals indicated significant differences between Lineages (Λ = 0.165, df error = 2847.8, P < 0.0001). No distinct clusters, however, were seen in canonical plots, which may be explained by the extensive overlap between most lineages in pairwise comparisons (Fig. 4).
Initially, a correct assignment of individuals to their lineage of origin was achieved in 72.8% of cases, but a considerably smaller classification success rate was obtained using leave-one-out cross-validation (LOOCV) Table 3. We observed similar patterns under different combinations of k and c between 0.1 and 0.5 (data not shown).
We conducted regularized DFAs (k = 0.1, c = 0.1) assigning individuals to the localities of origin for each lineage separately. All lineages and localities were used except Clade F, as it only consisted of one locality (Puerto Vallarta). All DFAs were significant (Table 4), but the percentage of individuals correctly assigned to their locality within each clade after LOOCV was low (range: 24.3-61.7%; Table 4), compared to before LOOCV (range: 97-100%).

Tests of phylogenetic signal
We did not detect widespread evidence of phylogenetic structure in the shape variables using either k or K statistics, but did notice a specific instance of structure related to the fifth shape PC (Table 5). This last result, however, appeared to be spurious, as no obvious differences were observed between lineages upon inspection. Although the statistical power of these tests is maximized when N > 20 (Pagel 1999;Blomberg et al. 2003), Pagel's k values are robust to the number of taxa included, whereas Blomberg's K values decrease as additional taxa are included (M€ unkem€ uller et al. 2012). Exploratory analyses incorporating guide trees with major clades subdivided into component lineages produced results consistent with these expectations and concordant with the results presented herein. We do not present these results, as the statistical significance of K values cannot be inferred using unresolved guide trees (Kembel et al. 2010).
K2P genetic distances ranged from 0.00 to 28.5% (mean = 20.2%, median = 21.6%). Euclidean distances in the morphological dataset ranged from 0.008 to 0.09 (l = 0.034, Median = 0.033). Lower dispersal of Euclidean distance values was observed at smaller genetic divergences. Regression of pairwise morphological distances against pairwise K2P genetic distances (Fig. 5) was significant (F = 79.0491, df error = 1375, P < 0.0001); however, the R 2 value suggested a poor fit between the data and the model (R 2 = 0.054). Mantel tests are known to be afflicted by high type-I error rates (Lapointe and Legendre 1995;Oberrath and Bohning-Gaese 2001;Nunn et al. 2006), and their use in phylogenetic comparative analyses has been discouraged (Harmon and Glor 2010). Nonetheless, we incorporated Mantel tests to determine whether the combination of all shape variables produced patterns different than those seen by evaluating shape variables independently.

Visualization of shape variation in Ligia
We present thin-plate-spline transformations for all major clades using vector magnifications of 109 for ease of comparison (Fig. 6). Because such magnifications have as much to do with statistical power as with the magnitude of effects, we selected individuals from each clade with the highest canonical score and provide those images as supplemental material (Fig. S1). Between-clade differences appear most pronounced in the cephalon, pleotelson, and midbody regions. Ligia individuals in Clade A exhibit an enlarged cephalon with small eyes. Their pleotelson is compressed, with the lateral posterior and the distal-most point almost parallel. Clade C has a pleotelson similar to that of Clade A; however, a somewhat rectangular body shape and small eyes with no enlarged cephalon may distinguish individuals from Clade C. Individuals in Clade B are characterized by an oval body shape with a normally sized cephalon and medium-sized eyes. As in Clade A, the pleotelson is compressed; however, the distal-most point protrudes well beyond the lateral posterior points. Clade D exhibits a slight invagination in the midbody region and medium-sized eyes on a regular cephalon. Their pleotelson appears less compressed than other clades, with the exception of Clade E. Although exhibiting a similar   pleotelson, Clade E has no midbody invagination. Also, individuals from this clade exhibit medium-sized eyes with a slightly larger cephalon than the rest of the body (i.e., the body tapers posteriorly). Clade S exhibits a body shape similar to those in Clade E; however, the body does not appear to taper, and the distal point of the pleotelson appears to protrude more extensively than in Clade E. Clade F specimens have very large eyes, with a drastic invagination in the segments prior to the pleotelson. Clade N has small eyes, with an oval body shape, a noncompressed pleotelson, and a large 1st segment.
We also present thin-plate-spline transformations for the largest and smallest individuals for both the overall dataset (i.e., Size effect) and for each lineage (i.e., Lineage*Size) (Fig. 7). In general, larger individuals exhibit a broader body and smaller eyes (relative to the total body size), with a distal point of the pleotelson that is slightly more protruding. All lineages appear to exhibit similar patterns, differing mostly in the magnitude of the effect. Individuals in clades A, C, D, and E exhibit the most obvious allometric effects (Fig. 7). Much subtler differences are observed between the large and small individuals in clades B, N, and S (Fig. 7).

Discussion
Taken together, our results suggest that L. occidentalis s. l. experiences morphological stasis. We detected significant differences in body shape among the major clades (variable Lineage) of L. occidentalis s. l., and among populations within these clades, but only the effect size of Lineage appears to be important (i.e., g 2 P > 0.2). Nonetheless, neither lineages nor populations can be discriminated on the basis of body shape, as correct classification rates of cross-validated DFAs were low (overall among lineages = 57.9%; range among populations within lineages = 24.3-61.3%). Such low rates of correct classification may be attributable to a large overlap of body shape among clades, as revealed by the canonical plots (Fig. 4). The much higher rates of correct classification prior to cross-validation implies that they were inflated, possibly reflecting over-fitting (Kovarovic et al. 2011). Failure to perform DFA cross-validation can produce severely biased results (Kovarovic et al. 2011) and thereby lead the investigator to erroneously conclude that the lineages under question can be distinguished with high accuracy. Unfortunately, studies that employ DFA often fail to report cross-validated results (e.g., Adams et al. 2014). Our results underscore the importance of cross-validation when conducting DFAs.
Extensive overlap among clades suggests body shape variation in L. occidentalis s. l. is highly constrained. It is possible that body shape variation can only occur within a narrow morphospace, henceforth the overlap among clades; and/or that morphological differentiation does occur, albeit at a very slow rate compared to genetic differentiation. Although the small slope of the Euclidean morphological distance versus genetic distance regression is suggestive of this pattern (Fig. 5), the poor fit of the model prevents robust inferences. Nonetheless, the largest Euclidean distances were observed at larger genetic distances. Due to the deep genetic divergences among the main lineages of L. occidentalis s. l., with some splits probably dating to the Miocene (Hurtado et al. 2010), the absence of strong body shape divergence is unlikely explained by a lack of genetic variation stemming from failure of mutations to arise. A similar finding of significantly different, but not fully diagnostic body shapes, was observed among three genetically divergent lineages of L. hawaiensis, an endemic of the Hawaiian archipelago . Although divergences among the three L. hawaiensis lineages are smaller than the deepest divergences among L. occidentalis s. l. lineages, they probably represent several million years of separation. Other cryptic species complexes of Ligia that include highly divergent lineages, which also appear to represent millions of years of separation, have been discovered in the Caribbean (Santamaria et al. 2014) and the Mediterranean (Hurtado et al. unpublished results). Despite a greater phylogenetic depth, considerable overlap in body shapes occurs between L. occidentalis s. l. versus L. hawaiensis (Fig. S2). Morphological stasis, thus, appears to be a common phenomenon within Ligia; a group that is at least $ 110 Maold based on the fossil record (Broly et al. 2013). Most Ligia species are coastal (Taiti et al. 2003;Broly et al. 2013) and exhibit morphological, physiological, and behavioral characteristics that are intermediate between ancestral marine and fully terrestrial isopods (Carefoot and Taylor 1995). Ligiidae occupy the most basal position within Oniscidea (Erhard 1998;Schmidt 2008), and it is believed that the ancestor of oniscideans had ligiid-like characters (Carefoot and Taylor 1995;Schmidt 2008). The traits that have allowed Ligia isopods to persist and diversify in their supralittoral environment may have contributed to constrain body shape evolution.
Strong stabilizing selection on body shape could be exerted by one or more inherent features of the intertidal habitat occupied by L. occidentalis s. l. (and other coastal Ligia). This selection must be efficient enough to counter the effect of drift, which is expected to be strong in many populations, due to their apparently small size and genetic isolation (Hurtado et al. 2010;Eberl et al. 2013). These isopods occupy a very narrow vertical range of rocky shores (i.e., from the water line to the supralittoral), where they track the shifting water line, and actively avoid the open sea. Because of their extremely low desiccation resistance, they must remain close to water, where they can take up water from droplets and puddles by capillarity and from water vapor directly from the air (Carefoot and Taylor 1995). They tend to be more active at night and remain hidden under rocks and in  crevices during much of the day (Carefoot and Taylor 1995;Hurtado et al. 2010), presumably to minimize water loss and avoid terrestrial/aerial predators (e.g., lizards, birds, and invertebrates; Hurtado et al. 2010;Grismer 1994;Brusca 1980). Their coloration enables them to blend in with the rocky substrate. Locomotion of these isopods is well adapted to the rocky substrate, where they dash for cover when threatened. On sand, however, their locomotion is extremely limited and their coloration is more conspicuous, rendering them highly vulnerable to predators. In addition, they may be more susceptible to wave action and desiccation, due to the lack of cover, in sandy substrate. Despite variation in the nature of the rocky beach type (i.e., gravel, cobbles, pebbles, boulders, and rocky bench; all of which are occupied by L. occidentalis s. l.), the rocky substrate may represent a rather homogeneous habitat that promotes the retention of a highly conserved body morphology. Accordingly, L. occidentalis s. l. appears to exhibit niche conservatism, at least in terms of substrate type (rocky) as it relates to body shape. Due to its effect on locomotory function, substrate type can be a critical determinant of morphology (e.g., Losos et al. 1997;Vervust et al. 2007;Goodman et al. 2008). Tylos, another cosmopolitan supralittoral endemic isopod genus, but that is restricted to sandy substrates, also exhibits high levels of cryptic diversity , and presumably body shape stasis. Therefore, this phenomenon may be common among permanent members of the supralittoral environment.
Body shape variation in terrestrial isopods (i.e., Oniscidea) is suggested to be constrained by a limited number of constructional pathways (Schmalfuss 1984). Accordingly, most of the~3700 species terrestrial isopods (Sfenthourakis and Taiti 2015) can be classified into five major functional categories of skeletal construction, which are correlated to ecological strategies and behavioral patterns: runners; clingers; rollers (e.g., Tylos); spiny forms; and creepers (Schmalfuss 1984). Runners dash when threatened to quickly hide under rocks or in crevices and are characterized by a narrow body, smooth tergites, long and strong pereopods, and a convex hind-margin of first epimeres. Tropical and subtropical species of Ligia, such as L. occidentalis s. l., are typically runners, whereas some temperate species of Ligia are suggested to be clingers (Schmalfuss 1984). Clingers, when threatened, tend to remain motionless clinging tightly to a solid substrate. Clingers differ from runners by a broader body, shorter pereopods, and a concave hind-margin of first epimeres (Schmalfuss 1984). Some body shape differences, thus, appear to occur between these two ecomorphs in Ligia.
Similarities in body shape among divergent lineages of L. occidentalis s. l. may be influenced by phylogenetic relatedness, geographic proximity (which may imply exposure to more similar environmental/ecological conditions), or stochasticity. We note, however, that phylogenetic relatedness and geographic proximity are highly confounded in L. occidentalis s. l., with phylogenetically closer lineages generally having adjacent distributions (Eberl et al. 2013;Hurtado et al. 2010;Hurtado et al., unpublished results). Clade A is sister to Clade BCDE; relationships within this last group are (B (E (C D))). Similarly, Clade N is sister to Clade S, and both occupy adjacent distributions in the Gulf of California, the former in the northern Gulf, whereas the latter in the southern Gulf. Clade F is a highly divergent lineage, whose distribution overlaps with that of Clade S. The relationships among the highly divergent ABCDE, NS, and F clades are uncertain.
Even though phylogenetic relatedness and geographic proximity are confounded, data on the clades to which misclassified individuals were assigned (Table 3) suggest that geographic proximity influences the degree of body shape similarity among clades. Individuals from geographically adjacent clades appear to have more similar body shapes. For example, the majority of the misclassified individuals by cross-validated DFA from Clade A were placed in Clade C, and vice versa. Most of the localities sampled from clades A and C are in the Northern Channel Islands. Therefore, geographic proximity appears to have a stronger influence than phylogenetic relatedness on the misclassification of individuals from clades A and C (A is sister to BCDE). Similar situations are observed for misclassified individuals of the other clades with the exception of Clade E. A noteworthy case is that of the misclassified individuals of Clade F, which were classified as members of Clade S in equal proportion to the correctly assigned individuals (i.e., 44.4%; Table 3). Clades F and S are phylogenetically distant, but the localities of Clade F examined are geographically nested within those of Clade S. The pairwise canonical plots show a remarkable separation between Clade F and all other clades, except Clade S, with which it exhibits complete overlap on the two major dimensions of among group variance (Fig. 4). Another example is that of Clade S, for which most misclassified individuals were assigned to geographically nearby clades N and E. Whereas clades N and S are sister lineages, Clade E is distantly related. Misclassification of individuals in Clade E, however, appears to be more stochastic, with few individuals being incorrectly assigned to the geographically adjacent Clade D. The lack of strong signal from phylogeny and genetic distance on body shape variation is consistent with the apparent effect of geographical proximity described above.
The apparent influence of geographic proximity on body shape similarity suggests, however, that the environment might impose at least some weak directional selection on shape variation. For example, the similarity between clades A and C, most of which were sampled on the Northern California Channel Islands, despite the marked differences in SST among insular localities occupied by the two clades (Eberl et al. 2013), suggests that insular ecological factors may be relevant to body shape (e.g., different or fewer terrestrial predators may be present in the islands). Determining the influence of extrinsic and intrinsic factors on body shape variation in L. occidentalis s. l. will require studies of ecological parameters, as well as of the genetic architecture of body shape, including an assessment of phenotypic plasticity (Wake 1991;Schlichting and Pigliucci 1998). This would allow evaluation of mechanisms other than stabilizing selection on body shape per se that may limit the evolution of this trait in L. occidentalis s. l. (e.g., genetic and developmental constraints, and stabilizing selection on traits correlated to body shape; Futuyma 2010).
We detected allometric effects on the overall body shape of L. occidentalis s. l., with g 2 P values (Table 2) suggesting body size to be the strongest determinant of overall body shape. In general, larger Ligia individuals (usually males) exhibit a disproportionately wider body than smaller individuals, a pattern reported for L. pallasii (Carefoot 1973). In L. hawaiensis, however, larger individuals exhibit a more elongated body with relatively smaller head than smaller individuals ). Thin-plate-spline visualizations suggest widening of the body is not uniform across L. occidentalis s. l. lineages. Developmental (Stern and Emlen 1999) or ecological differences (Pfennig 1992) may be responsible for the differences in the magnitude of this effect observed in L. occidentalis s. l. Three of the lineages (C, D, E) exhibiting the deepest widening of the body constitute a wellsupported monophyletic clade, with those exhibiting no obvious allometric effects (N and S) also forming a monophyletic group (Hurtado et al. 2010). These patterns may be indicative of the shared evolutionary history of these lineages and may be due to shared developmental constrains. On the other hand, environmental factors may be at play. Growth rates of isopods are known to be affected by environmental factors such as temperature (Strong and Daborn 1980;Holdich and Tolba 1981;Donker et al. 1998), food availability (Reichle 1968), and exposure to pollutants (Donker et al. 1993). In general, lineages in the colder Pacific Ocean (A, C, D, E) exhibited obvious allometric effects, whereas those in the warmer Gulf of California (N, S) did not. Furthermore, we observed some obvious allometric effects in some Clade N localities (N7, N9, N12), suggesting that differences in allometric effects on body shape may also correspond to ecological factors and not phylogenetic trajectory. As the distribution of L. occidentalis s. l. lineages closely matches changes in sea surface temperatures (Eberl et al. 2013), additional work is needed to establish the contributions of ecological differences and phylogenetic relatedness on the observed differences in the magnitude of allometric changes.
Ligia occidentalis s. l. appears to truly represent a hypercryptic species complex, lacking diagnostic morphological differences among putative species. Preliminary examination for differences among lineages based on traditional characters used in Ligia has failed to reveal obvious differences. In addition, one of us (CAS) has preliminarily examined SEM micrographs of the appendix masculina and 7th pereopod of individuals representing the main clades, and has not found clear differences. The apparent hypercryptic nature of this taxon brings about challenges for the conservation of a unique and highly vulnerable biodiversity. The current recognition of a single widely distributed species hampers preservation efforts. Unique cryptic lineages have very restricted geographic ranges, usually constrained to discrete beaches, with some that appear to have small population sizes. The rocky supralittoral habitat in the range of L. occidentalis s. l. is very vulnerable to destruction and modification by human activities (e.g., construction and pollution), which have increased as human populations and tourism expand (Hurtado et al. 2010; and populations of this isopod across its range have been shown to accumulate toxic contaminants from anthropogenic activities (Garc ıa-Hern andez et al. 2015). Storms and hurricanes may sweep local populations (Hurtado, personal observation), whereas, in the long term, sea level rise could also threaten the persistence of populations. In the absence of morphological differences that can be taxonomically diagnostic, we suggest the use of molecular information to taxonomically classify L. occidentalis s. l. in a way that reflects and helps protect its unique diversity. Defining species based on molecular data alone has been conducted in other cryptic species complexes (Dumas et al. 2015;Wade et al. 2015). In L. occidentalis s. l., the process may be facilitated by the discrete geographic distribution of lineages. We suggest the use of a divergence cutoff to assign species and subspecies, which should aid in the conservation of the rich biodiversity found within this clade.
Dan Richards (Channel Islands National Park), Channel Islands National Marine Sanctuary, Darsee Guttilla (The Catalina Conservancy), Mary Wicksten (Texas A&M University), Renate Eberl, and La Cooperativa de Pescadores Nacionales de Abul on in Isla Cedros. Ellen Giddens, Maureen Frank, and Micah Wagner helped with laboratory work.
This work was supported by National Science Foundation (NSF) grant DEB 0743782 to LAH and MM, by a Texas A&M University-Consejo Nacional de Ciencia y Tecnologia, Mexico (CONACyT) grant to LAH, and by Texas A&M University. Partial support for CAS was provided by the Sloan Minority PhD Program and by the Hispanic Leaders in Agriculture and the Environment Program, and by the Tom Slick Senior Graduate Fellowship. The open access publishing fees for this article have been covered by the Texas A&M University Online Access to Knowledge (OAK) Fund, supported by the University Libraries and the Office of the Vice President for Research.