The First Insight into the Patterns of Size and Shape Variation of a Microcerberid Isopod

: Cryptic species are a biological phenomenon only recently recognized due to progress in molecular studies. They pose a signiﬁcant challenge to conventional taxonomic work since these species manifest low morphological differences, but considerable genetic disparity. New taxonomic methods are in development but have yet to be tested for many animal groups. Isopods belonging to the suborder Microcerberidea are one such group. The Asian microcerberid isopod, Coxicerberus fukudai (Ito, 1974), is a major component of marine interstitial fauna with suspected cryptic species inhabiting Japan and Korea. We chose six Korean populations with high molecular interpopulations divergence and applied 2D landmark-based geometric morphometrics to cephalic sensilla, pleonal points, and male pleopod II. This quantitative approach allowed us to study interpopulation size and shape variations, morphospace structure, and whether the morphological pattern mirrored the genetic species. We determined that a high degree of interpopulation size variation signiﬁcantly inﬂuences shape changes. Once we removed the allometric effect, the size-corrected male pleopod II shape variations yielded a new species, C. jangsaensis sp. nov. At the same time, we were able to resolve the C. fukadai species complex.


Introduction
The body shape of a species often reflects the accumulation of historical contingencies naturally selected by the environment [1]. Certain environmental conditions can impose morphological stability resulting in population lineages, living in allopatry, with no apparent morphological differences [2]. This phenomenon, called "cryptic species", refers to a group of closely related genetic species morphologically indistinguishable and erroneously classified as a single nominal species [3,4]. This unexpected genetic diversity has greatly altered not only global species diversity estimates, but conventional taxonomy as well [5,6]. An integrative approach combining morphology and molecular data might resolve these taxonomic problems [7][8][9][10][11][12][13]. This is because applying multiple lines of evidence is more likely to better define species boundaries [8]. Quantitative analysis is a complementary approach for providing objective data along with traditional taxonomic work. Landmark-based geometric morphometrics (LBGM) and an outline contour, for instance, can analyze complex biological structures using landmark sets, which are located on every specimen [11]. This technique analyzes shape data defined by a superimposed landmark configuration with size, rotation, and translation effects removed [14]. Shape variations are presented as statistical components describing landmark shifts. These can be visualized using a wire frame [15]. Landmark-based geometric morphometrics can detect subtle variations, providing insight into the asymmetry [16][17][18][19], intrapopulation shape variations [20][21][22], as well as ontogenetic variation [23,24]. Applying LBGM techniques with respect to cryptic species is an emerging taxonomic field, especially in microcrustaceans. Several morphological structures have been used for this purpose: carapace outline in ostracods [25,26], cuticular pores in copepods [27,28], epimera in amphipods [29], and body surface in cladocerans [30]. Only a few terrestrial oniscidean species have been studied in the order Isopoda, and those were to discern morphological and genetic diversity. Results showed that allopatric populations of the common sea slater, Ligia Fabricius, 1798, display a high level of genetic divergence but no morphological differences between coastal populations [31,32]. Subsequent geometric morphometrics studies detected significant morphological variation between geographical lineages in lateral tergite tips of each body segment [33,34].
Members of the isopod suborder Microcerberidea are small animals (less than 1 mm) inhabiting marine and freshwater interstitial habitats [35]. Since Karaman (1933) [36] described the first microcerberid species, nearly 30 species have been recognized so far [37]. The group is morphologically homogenous, lacking pigmentation and having an elongated body, adapted to interstitial water environmental conditions [38]. Despite recent phylogenetic and ecological studies [39,40], there are many biological aspects that remain largely unknown. For example, intraspecific morphological variation studies are almost absent even though some microcerberid species are widely distributed geographically [41][42][43][44]. Most taxonomic studies use a small number of individuals, which makes variability detection difficult. Kim et al. (2018) [40] studied variability within the East Asian species, Coxicerberus fukudai (Ito, 1974), originally described from northern Hokkaido's sand beach and later found in Korea. Molecular analyses revealed two divergent Korean lineages, one distributed in the northern part of South Korea, and forming a monophyletic clade with C. fukadai, and the other, separate lineage, to the southern part. Kim et al. (2018) [40] observed morphological abnormalities in antenna among southern population individuals, but detected no consistent morphological disparity between the two lineages. A high genetic divergence (mtCOI: 10.5%; Cytb: 9.5%) and absence of consistent morphological differences indicates the presence of cryptic species. Male genitalia morphology in microcerberid isopods is a major diagnostic characteristic [38,44], so the lack of a difference between genetically distinct populations is unusual and may indicate a new impediment to traditional microcerberid taxonomy.
In this study, we used LBGM to evaluate morphological variability between two genetically divergent Korean populations. We chose three structures: cuticular sensilla on the dorsal surface of the cephalothorax, the pleon's anatomical points along each pleonite, and the male's pleopod II. Sensilla [45,46] refer to minute cuticular sensory receptor organs covering a vast area of the crustacean body [47]. Due to their variability in size, shape, number, and position, many studies have investigated their taxonomic importance in various crustacean taxa [48][49][50][51][52][53][54]. The pleon's anatomical points along each pleonite and pleotelson boundary are used to estimate segment length ratio changes, which is considered important in the taxonomy of isopods [33,34]. The male pleopod II endopodite displays a particularly modified morphology, so-called "arm and hammer". This is a synapomorphy of Asellota and Microcerberidea, phylogenetically closely related groups, and supposedly it has a constrained evolution due to interlinked parts and its role as a continuous sperm conduit [55]. The endopodite of the male's pleopod II is specialized for conveying sperm packets to female genitalia and has a variety of forms that have been generally sufficient for species identification. Our initial study could not detect differences in the anatomy of male's pleopod II between the two lineages of the Korean populations of C. fukudai. Nevertheless, we chose this structure for our quantitative analysis based on the fact that mismatch in the genitalia has long been posited as a mechanism of reproductive isolation [56]. The primary aim of our study is to resolve the C. fukadai species complex with LBGM, as a novel approach in the microcerberid and isopod studies in general.

Sample Collection
We collected samples from Gazin, Sokcho, Dongho (lineage A), Detan, Jangsa, and Imrang (lineage B) beaches in Korea for our two genetic lineages comparative analysis study [40]. We collected samples on two occasions in September and October 2020 by filtering the water with a hand net (mesh size 46 µm) and immediately preserved the samples in 99% ethanol. We kept all collected animals in three plastic bottles (200 mL) separately from each beach and preserved them in the portable refrigerator until reaching the laboratory. All newly collected animals with geographical information were deposited in the collection of the National Institute of Biological Resources (NIBR), Incheon, Korea (Table S1).

Genetic and SEM Analyses
To confirm the population's previously established genetic identity, we randomly selected 1-3 individuals from each sample and extracted their DNA. We amplified the mitochondrial cytochrome b sequences following the same procedure as Kim et al. (2018) [40]. We confirmed all obtained sequences with BLAST search [57], which were deposited in Gen-Bank (Table S2) [58]. We performed sequence alignment using the MAFFT v7.313 [59] with those of C. fukudai already published and publicly available. We used the TCS algorithm implemented in PopART [60] to visualize the haplotype network for the aligned sequence dataset. We evaluated the best substitution model (TIM1) as well as maximum likelihood (ML) trees using the iqtree version 1.6.9 [61] with 1000 ultrafast bootstrap replicates. We constructed an ML phylogram with iqtree default parameters. We transferred specimens for the observation under the scanning electron microscope (SEM) to isoamyl acetate for 20 min and dried them in a critical-point dryer (Hitachi E-1010; Hitachi, Tokyo, Japan). We mounted dried specimens on an SEM stub coated with gold using a sputter coater to a thickness of 15-30 nm. We then photographed the coated specimens with a Hitachi S-3400 scanning electron microscope at Chungang University, Seoul, Korea.

Morphological Data Acquisition
We confirmed individual maturity using a combination of the following characteristics: (1) length of the pereonite VII (longer in adults vs. much shorter in juveniles); (2) length of the pereopod VII (incompletely developed in juveniles); (3) length of the male pleopod II (absent in juveniles); (4) whole body length (adults generally exceed 0.8 mm). We confirmed individual sex using a combination of the following characteristics: (1) pereonite VII penial papillae (present in males vs. absents in female); (2) male pleopod II (present in males vs. absent in females); (3) whole body length (males measure less than 1 mm vs. females longer than 1 mm). Information on individual animals is summarized in Table S2. We dissected the animals using a stereo dissecting microscope, Olympus BX51 (Olympus, Tokyo, Japan). We separated the cephalothorax and pleon from the body and mounted them immediately on microscope slides with lactophenol and two coverslips. We made line drawings of cephalothorax, pleon, and male's pleopod II endopodite using a compound microscope, Olympus BX51 (Olympus, Tokyo, Japan), equipped with a camera lucida. To make a dorsal perspective angle perpendicular to the microscope objective, we manipulated the coverslip by hand, and focused on the symmetrical view [28]. After each specimen was adequately positioned consistently, we made pencil drawings of the cephalothorax, pleon, and pleopod II endopodite. We scanned the images and generated a TPS file using tpsUtil software (File S1) [62]. We digitized the chosen landmark (LM), all Type I [63], twice using tpsDig2 software [64] to estimate digitization-related errors [65]. We analyzed the sensilla distribution (SD) on cephalothorax, pleonal points (PP), and pleopod II endopodite (PE) as separate datasets. We selected eight sensilla ( Figure 1A): two located above a pair of incisions on the dorsoposterior margin (LMs 1 and 2), two along the posterolateral margins (LMs 3 and 4), two on the medial margin (LMs 5 and 6), and two along the anterolateral margin (LMs 7 and 8). We selected eight anatomical points of pleon ( Figure 1B): two on the right and left vertex of the pleonite I anteriorly (LMs 1 and 2), two on the right and left vertex of the pleonite II anteriorly (LMs 3 and 4), two on the right and left vertex of the pleotelson anteriorly (LMs 5 and 6), and two points on the pleotelson's posterior end (LMs 7 and 8). The endopodite is an asymmetric object with overall shape out curved, its distal end is bifurcated, of which the lateral projection is elongate and distally pointed. We noticed from the preliminary observation that the location of lateral projection seems to vary across the populations, which was previously unreported. We, therefore, selected six anatomical points at the distal end of right PE at the ventral position ( Figure 1C,D) to evaluate the relative position of lateral projection against tip of medial body: one on the distal tip of medial body (LM 1), one on the proximal groove's vertex between the body and the projection (LM 2), one on the medial side of the projection symmetrically in the same position as the distal tip of the body (LM 3), one on the distal tip of projection (LM 4), one on the lateral side of the projection symmetrically in the same position as the distal tip of the body (LM 5), and one on the proximal end of the projection symmetrically in the same position as the LM 2 (LM 5). The SD, PP, and PE datasets included 360, 370, and 182 digitized images (from the original 180, 185, and 91 pencil drawings), respectively. We summarized the image numbers generated according to the classifier (lineage and population, and sex) in Table 1. Except when performing the Procrustes analysis of variance (ANOVA) to test measurement error, we used only the first digitization for the subsequent LBGM analyses.
Water 2021, 13, x FOR PEER REVIEW 4 of 18 two along the anterolateral margin (LMs 7 and 8). We selected eight anatomical points of pleon ( Figure 1B): two on the right and left vertex of the pleonite I anteriorly (LMs 1 and 2), two on the right and left vertex of the pleonite II anteriorly (LMs 3 and 4), two on the right and left vertex of the pleotelson anteriorly (LMs 5 and 6), and two points on the pleotelson's posterior end (LMs 7 and 8). The endopodite is an asymmetric object with overall shape out curved, its distal end is bifurcated, of which the lateral projection is elongate and distally pointed. We noticed from the preliminary observation that the location of lateral projection seems to vary across the populations, which was previously unreported. We, therefore, selected six anatomical points at the distal end of right PE at the ventral position ( Figure Table 1. Except when performing the Procrustes analysis of variance (ANOVA) to test measurement error, we used only the first digitization for the subsequent LBGM analyses.

Geometric Morphometric Analyses
We used algorithms implemented in Morpho J package software version 1.07a [66] for all LBGM analyses. We aligned and superimposed all landmark configurations with generalized Procrustes analysis (GPA) to remove the effects of non-shape variation [14]. Of SD, PP, and PE datasets, the former two ones were object symmetry [65], with the axis of plane passing through the LM configuration; each LM has a corresponding pair on the opposite side ( Figure 1A,B). The GPA also detects directional and fluctuating asymmetry and shape differences. For this study, however, we considered the symmetric component only, which represents average shape changes between the left and right sides. Because the potential discrepancies in angle can be induced by specimen positioning on the microscope that has significant effects on the asymmetric component [28,67], we converted the Procrustes shape coordinates into a covariance matrix [68]. As a size proxy, we estimated the centroid size (CS) for each individual from the raw LM coordinates (Table S3) [15]. We calculated the CS as the square root of the sum of squared distances for a set of centroid LMs [69]. Descriptive statistical values of CS such as means, standard deviation, and standard errors were calculated and displayed with box plots generated by a web-tool for generation of box plots, "BoxPlotR" [70].
We analyzed the SD, PP, and PE covariance matrices independently (File S2). We used Procrustes ANOVA to test group structuring evidence in the overall dataset with lineage and sex as classifiers as well as the digitizing errors [71]. We performed regressions of shape onto size to test allometry using regression scores and CS. [72,73]. The null hypothesis states that shape develops isometrically; thus, a statistically significant result demonstrates that shape changes with increasing size according to a predictable model [74]. We applied a permutation test [75] to assess the statistical significance against the null hypothesis. The number of shape variations determined by the regression was expressed as a percentage of the total variation percentage around the mean [24]. We calculated the SD, PP, and PE datasets residual components to subtract the portion of shape variation predicted by the regression for further analyses. We analyzed the residual shape component using principal component analysis (PCA), which is the most efficient method for examining the variation of multiple variables within a single sample. This analysis is often used as the first exploratory analysis of a large dataset composed of several samples to provide a visual impression of overall shape variation. [69]. We performed the PCA twice on the entire dataset and then on datasets of females and males separately with the wire frame to visualize SD, PP, and PE's average shape variation along major PCA axes. We used the LM displacement described in the PCA to establish the new species' taxonomic diagnosis. We employed PE dataset residual components for discriminant function analysis (DFA) and canonical variate analysis (CVA). Both are multivariate LBGM methods that produce a criterion for reliably distinguishing two lineages and among multiple populations, respectively. DFA and CVA results indicate the separation among groups by maximizing distances between group means relative to the variation within groups [76]. For the result of DFA, we applied cross-validation to test the reliability of groups separation [76]. We used multivariate statistics as Mahalanobis distances [77] to estimate the difference between groups. The permutation test assessed the statistical significance against the equal group means' null hypothesis. We finally used the cross-validation score of DFA and the Mahalanobis distances of CVA to evaluate the availability of pleopod II endopodite for delimiting the groups.

Confirmation of Genetic Indentity
The haplotype network ( Figure 2A) detected six haplotypes, with the Gazin, Sokcho, and Dongho populations (lineage A) sharing a haplotype separated from the lineage B's ones with over 30 mutational substitutions. The Imarang population showed four haplotypes; each of them is separated by one mutational substitution. The Imrang one was the only location where multiple haplotypes exist. The Detan, and Jangsa, populations shared a haplotype distinguished by one substitution from the Imrang's one. The ML phylogram ( Figure 2B) revealed two distinct clades with absolute certainty. The Gazin, Sokcho, and Dongho sequences formed a distinct clade representing a diverged lineage from that of Detan, Jangsa, and Imrang's clades. All Imrang sequences formed separated clades from the Detan and Jangsa ones. to provide a visual impression of overall shape variation. [69]. We performed the PCA twice on the entire dataset and then on datasets of females and males separately with the wire frame to visualize SD, PP, and PE's average shape variation along major PCA axes. We used the LM displacement described in the PCA to establish the new species' taxonomic diagnosis. We employed PE dataset residual components for discriminant function analysis (DFA) and canonical variate analysis (CVA). Both are multivariate LBGM methods that produce a criterion for reliably distinguishing two lineages and among multiple populations, respectively. DFA and CVA results indicate the separation among groups by maximizing distances between group means relative to the variation within groups [76].
For the result of DFA, we applied cross-validation to test the reliability of groups separation [76]. We used multivariate statistics as Mahalanobis distances [77] to estimate the difference between groups. The permutation test assessed the statistical significance against the equal group means' null hypothesis. We finally used the cross-validation score of DFA and the Mahalanobis distances of CVA to evaluate the availability of pleopod II endopodite for delimiting the groups.

Confirmation of Genetic Indentity
The haplotype network (Figure 2A) detected six haplotypes, with the Gazin, Sokcho, and Dongho populations (lineage A) sharing a haplotype separated from the lineage B's ones with over 30 mutational substitutions. The Imarang population showed four haplotypes; each of them is separated by one mutational substitution. The Imrang one was the only location where multiple haplotypes exist. The Detan, and Jangsa, populations shared a haplotype distinguished by one substitution from the Imrang's one. The ML phylogram ( Figure 2B) revealed two distinct clades with absolute certainty. The Gazin, Sokcho, and Dongho sequences formed a distinct clade representing a diverged lineage from that of Detan, Jangsa, and Imrang's clades. All Imrang sequences formed separated clades from the Detan and Jangsa ones.

Testing Variations in Size and Shape
ANOVA results (Table 2) yielded a negligible digitizing error for all datasets with the individual variability mean squares (MS) and F values by far exceeding the error values. Lineage was a significant contributor to overall shape and size variations for SD and PE datasets. It was highest in the PE shape component dataset (F = 29.16, p < 0.0001), followed by the PE size dataset (F = 13.48, p < 0.0001), SD size (F = 11.96, p = 0.0007), and shape dataset (F = 5.18, p < 0.0001). The lineage did not contribute to the size and shape in the PP datasets (size F = 0.65, p = 0.4217; shape F = 1.5, p = 0.1746), while the sex contributed significantly to the PP dataset (size F = 89.17, p < 0.0001; shape F = 54.2, p < 0.0001), and the SD shape dataset (F = 4.03, p = 0.0005). In the SD dataset, sex has not contributed to the size component of variability. The box plots summarize differences in the CS in detail and emphasize substantial sexual size dimorphism. In the SD dataset ( Figure 3A), the largest individual was in the Gazin females and the smallest one was in the Jangsa females. The Gazin, Sokcho, and Detan individuals generally showed larger CS indicated by interquartile ranges than the Jangsa, Dongho, and Imrang ones. Sexual dimorphism was, however, not much comparable as most populations displayed similar interquartile ranges between sexes, with some exceptions from the Gazin and Dongho showing larger female CS. On the other hand, the PP dataset ( Figure 3B) showed the Gazin females and Dongho males comprising the largest and smallest individuals. The sexual dimorphism was a general trend in all populations, but it was more pronounced in the Dongho, Gazin, and Jangsa populations than others in terms of the absence of the overlapped ranges. PE dataset (not presented) displayed the CS variation ranging from four to nine across the populations. Sokcho and Imrang individuals included the largest and smallest individuals, respectively. Table 2. Size and shape variation of sensilla distribution, pleonal points, and pleopod II endopodite inferred by Procrustes ANOVA where a randomized permutation procedure was applied (10,000 iterations). F: Goodall's F critical value, p: probability of finding a randm value larger than the observed value.

Allometry and Size-Corrected Shape Variation
The SD entire dataset regression for isometric shape development estimated a significant allometric effect (16.76%, p < 0.0001), so the null hypothesis was rejected. The PCA based on the residuals described the major shape variations of overall dataset with the first two axes carrying 57.4% (PC1 = 36.8%; PC2 = 20.6%) of the total variance. Despite a high individual variability, the PCA did not describe any difference between lineages. The sex-specific PCA using the residuals explained the major male (not presented) and female ( Figure 4A

Allometry and Size-Corrected Shape Variation
The SD entire dataset regression for isometric shape development estimated a significant allometric effect (16.76%, p < 0.0001), so the null hypothesis was rejected. The PCA based on the residuals described the major shape variations of overall dataset with the first two axes carrying 57.4% (PC1 = 36.8%; PC2 = 20.6%) of the total variance. Despite a high individual variability, the PCA did not describe any difference between lineages. The sex-specific PCA using the residuals explained the major male (not presented) and female ( Figure 4A  The PCA of female, however, did not show any difference in the shape between lineages ( Figure 4A) and between populations ( Figure 4B). The wire frame demonstrated that some Sokcho, Dongho, and Detan individuals showed the wider and shorter SD than average due to LM1 and LM2 shifted antero-medially, LM3 and LM4 shifted antero-laterally, and LM5, LM6, LM7, and LM8 shifted postero-laterally ( Figure 4C upper panel). About three Imrang individuals displayed the SD narrower and longer than mean shape due to LM1 and LM2 shifted postero-laterally, LM3 and LM4 shifted postero-medially, and LM5, LM6, LM7, and LM8 shifted antero-medially ( Figure 4C lower panel).   The regression of the PP overall dataset estimated a significant allometric effect (10.63%, p < 0.0001), rejecting the null hypothesis of isometric shape variation. The first two PCA axes based on the residuals accounted for 81.7% (PC1 = 70.2%; PC2 = 11.5%) of the total shape variation. However, lineages did not show any meaningful clustering in the morphospace. The first two PCA axes of the male (not presented) and female residuals analyses ( Figure 4D,E) carried 81.91% (PC1 = 68.07%; PC2 = 13.85%) and 79.02% (PC1 = 67.06%; PC2 = 11.97%) of the total variability, respectively. The PCA of both sexes showed hardly any distinguishable clustering among populations. Nevertheless, the PCA of the female dataset showed the morphospace clustering of Jangsa, Dongho, and Detan populations on one side (with the center of gravity in the positive part of the PC1) and the one of Gazin, Sokcho, and Imrang on the other (with the center of gravity in the negative part of the PC1). The wire frame of PC1 ( Figure 4F, upper panel) showed a narrow and elongated PP shape for the Jangsa, Dongho, and Detan populations due to LM1 and LM2 shifted antero-medially, LM3, LM4, LM5, and LM6 shifted postero-medially, and LM7 and LM8

Lineage Differentiation and Morphological Distances
shifted antero-medially. The wire frame of PC1 for the Gazin, Sokcho, and Imrang ones ( Figure 4F lower panel) differed from the former group by having LM3 and LM4 shifted antero-medially, and LM7 and LM8 postero-medially.
The PE dataset regression analysis estimated a significant allometric effect (6.02%, p = 0.002) so the null hypothesis of isometric shape development was also rejected. The PCA based on the residuals from the regression analysis showed the major shape changes occurring in the first two PCs, accounting for 83.7% (PC1 = 49.1%; PC2 = 34.6%) of the total variance. The PC1 emphasized morphological differences between lineages ( Figure 5A), with the lineage A having the center of gravity on the negative part of the axis, and lineage B on the positive. When populations were used as determinant of the PC1 morphospace clustering ( Figure 5B), Gazin, Sokcho, and Dongho individuals mostly occupied the negative space between −0.4 and 0. The Jangsa and Imrang individuals had variations ranging between −0.1 and 0.3 clearly with a positive center of gravity, while the Detan individuals had a relatively wider variation, ranging between −0.3 and 0.2 with a center of gravity near 0. The PC1 described the wire frame for the Gazin, Sokcho, and Dongho individuals ( Figure 5C, right panel), which differed from that of the latter group ( Figure 5C, left panel) by the following LMs displacement: LM1 and LM3 shifted proximo-medially, LM2 and LM6 shifted disto-medially, LM4 shifted disto-laterally, and LM5 shifted proximo-laterally.

Lineage Differentiation and Morphological Distances
The permutation test of the cross-validation of DFA and CVA rejected the PE null hypothesis for equal group means between lineages (p < 0.0001) and between populations (p < 0.0001). The Mahalanobis distances between lineage A and B was 2.2277 (p < 0.0001). The cross-validation of DFA based on residuals (Table 3)

Lineage Differentiation and Morphological Distances
The permutation test of the cross-validation of DFA and CVA rejected the PE null hypothesis for equal group means between lineages (p < 0.0001) and between populations (p < 0.0001). The Mahalanobis distances between lineage A and B was 2.2277 (p < 0.0001). The cross-validation of DFA based on residuals (Table 3) (Table 4) revealed that the Jangsa population was the most distantly related to lineage A, followed by the Imrang and the Detan ones. Lineage B showed a high intraspecific variation-the Jangsa individuals were accurately distinguished from the Detan ones (2.38, p < 0.0001), while those from Imrang were barely separated from the Jangsa (3.17, p = 0.005) and Detan (2.66, p = 0.04) ones. Etymology. The new species is named after the type locality. Diagnosis. General characteristics including habitus, appendages segmentation, and appendage setal formula as given by Kim et al. (2018) [40] in redescription of C. fukudai. Jangsa and Detan individuals, with pleopod II endopodite PC1 shape variation based on the following LM displacements: LM1 displaced disto-medially, LM2 proximo-laterally, LM3 disto-medially, LM4 proximo-medially, LM5 disto-medially, and LM6 proximo-laterally.
Remarks. One of the major differences between C.fukudai and C.jangsaensis sp. nov. includes the relative distance between LM1, LM3, LM4, and LM5. The Jangsa, the Imrang, and over half of the Detan individuals representing C. jangsaensis sp. nov. displayed a shape change with decreased distance between those landmarks, resulting in the proximity of the projection's distal tip and the medial body's one (see Figure 5C Right). The Dongho, the Gazin, and the Sokcho individuals of C. fukudai showed the shape change with increased distance between the landmarks, representing the remoteness of the projection's distal tip and the medial body's one ( Figure 5C, left). In addition, the relative movement of LM1, LM2, LM3, LM5, and LM6 described another distinct shape change between two species. The lengthen distance among those landmarks in C. jangsaensis sp. nov. showed overall a deep proximal rut represented by LM1-3, while those of C. fukudai showed the shortest distance, forming a relatively shallow groove. As a result, overall LM displacements indicates that C. jangsaensis sp. nov. tends to display the PE's lateral projection less developed distally than C. fukudai's one.

Discussion
One of the most significant results of the LBGM analysis is a high degree of CS differences in SD, PP, and PE datasets. The grouping factors play a dominant role in the overall size variations. In the SD and PE datasets, lineages A and B differed significantly from each other in size, while the PP dataset showed hardly any statistical differences. Sex was a significant factor accounting for the overall size variation in the PP dataset, but not in SD. Such size variation consequently yielded a significant allometric association, which is the main cause of the disparity or similarity between lineages and sexes. This is particularly obvious when comparing the shape variation in SD and PP datasets with and without allometry. The shape variation produced by a continuation of allometric association is not an independent morphological criterion and is invalid for the species delineation [78], but our discovery provides the first insight into the size variation in the microcerberid isopods. Initially, Kim et al. (2018) [40] reported sexual dimorphism in body size in all studied Korean populations of C. fukudai: the overall body length of adult females lies between 1 and 1.5 mm and is less than 1 mm for males, which is a common trend in other congeners [42,45,[79][80][81][82]. Given such interspecific variation in size, overall body length and sexual size dimorphism discrepancies in the sensilla distribution are prominent features of lineage A and lineage B. Secondly, the genetic lineages did not correspond to the similarity in CS given that all populations showed a significant difference in the intrapopulation variation. Since size difference can either be a result of sympatric speciation [83] or a temporary response to the ecological constraints [84][85][86], we need to explore the cause further. Potential associations between size variation and seasonal changes in temperature, nutrient levels, pH, salinity, and predation [87] would benefit from sampling over an extended period of time; microcerberid ecology is poorly known other than salinity and temperature tolerance ranges for one species [39].
The PCA analysis of regression residuals in the sensilla dataset revealed that lineages A and B have a similar intrapopulational variability ranges, resulting in indistinguishable morphospace clustering. This is interesting because results from other microcrustacean studies indicate that cuticular organs are important in taxonomy [27,53,[88][89][90][91]. Considering a deep genetic divergence between two lineages with a split likely occurring between 3.6 and 2.1 million years ago [40], the absence of shape differences among the populations cannot be explained by a lack of genetic variation. Instead, it is the result of stabilizing selections of cuticular pore distribution patterns. Selection of this sort is known to constrain populations to relatively constant and high-fitness phenotypes and thus limits the frequency of directional or disruptive evolutionary changes [92]. The genus Coxicerberus is a marine interstitial group, displaying morphological characteristics associated with the life in narrow, dark spaces [38]. Eyes are absent, forcing them to rely strongly on other sensory organs [93]. The cephalothorax's dorsal surface sensilla are assumed to play a prominent role in transmitting mechanical information to the central nervous system for efficient food sourcing and predator avoidance [94][95][96]. Despite their functional importance, we do not know to what extent these sensilla are species-specific since only a few Coxicerberus studies mention the number of sensilla along the cephalothorax's lateral margin [42][43][44]82,93]. More sensilla studies are necessary in order to more accurately address their evolutionary stability and interspecific differences. A recent LBGM copepod study describes how cuticular pores and sensilla distribution patterns can be used to distinguish between closely related species [97].
The PCA analysis of the pleon dataset regression residuals showed a significant individual variability in the pleonite I, II, and pleotelson, altering their L/W ratio. However, there is no morphological evidence to distinguish two lineages in both sexes. Instead, we detected a distinct population clustering, but only in the female dataset. Accordingly, Gazin and Sokcho populations clustered together (lineage A) and are separated from the Jangsa (lineage B), the former groups showing much broader and shorten pleons than the latter one. The majority of Detan clustered with the Jangsa population with high similarity in the pleon shape, corresponding to the genetic membership. The Imrang and Dongho individuals did not follow this pattern and displayed a broad and elongated pleon, respectively. This phenomenon may indicate that the female pleon underwent a local adaptation across the populations in a different direction from the male. The female-specific variation is worth studying further, as the microcerberid taxonomy has mostly relied on male features. Female morphology is generally homogenous across taxa, with the overall shape lacking distinct qualitative traits, consequently being neglected in previous studies. Given that the geometric trends can differ between sexes, it is essential to test a strong sexual dimorphism with molecular markers in order to confirm con-specificity. Marine isopod populations show a strong female-biased sex ratio [98][99][100], from just over 50% to 14:1 in some Jaera Leach, 1814 species [101]. We also observed a similar phenomenon for C. fukudai complex, which was problematic when comparing male specimens as they are morphologically indistinguishable or are not easily available [40]. Therefore, extrapolating our results regarding the importance of female characteristics to other taxa would be beneficial for more reliable species delineation. PCA, the cross-validation, and CVA analyses of the pleopod II dataset residuals showed a distinct separation of lineages A and B leading to a new species description. Coxicerberus jangsaensis sp. nov. displayed the shape variation of pleopod II endopodite's distal tip, which is apparently opposite to that of C. fukudai. Such variation in the genital parts is occasionally associated with the copulation success rate, resulting in the reproductive isolation among morphologically diverged individuals [56,102]. Wilson (1991) [55] described that the copulation of Asellota is mainly controlled by exopod holding the pleopod II and forcing it into the female's genitalia when sperm packets are conveyed. Coxicerberus jangsaensis sp. nov. shows the lateral projection of endopodite is structurally flexible and located near the opening where the appendix masculina is developed. Given the phylogenetic closeness and synapomorphy of the pleopod morphology between Asellota and Microcerberidea [35], we assume this body part interacts with the proximal part of the appendix masculina during the insemination. Indeed, the male pleopod II endopodite's distal tip is one of the most variable characters across the species that has been used to determine a morphological boundary or phylogenetic relatedness. According to the species grouping by Baldari and Argano (1984) [44], C. fukudai has a strong morphological affinity with C. abbotti Lang, 1961 and C. kiiensis Nunomura, 1973 in terms of sharing the overall shape of pleopod I endopodite, which is out curved, covered by rows of setule medially, and tapering distally. One of the most notable differences between those species is the morphology of the endopodite's distal tip, which is rounded/pointed without any projection in the latter two species. Based on this observation, we speculate that the geometric variability in the PE's distal tip is a significant morphological evidence accounting for the divergence and reproductive isolation between C. fukudai and C. jangsaensis sp. nov.
With a reverse taxonomy concept [103], we employed a quantitative LBGM approach using three morphological characters, pleon, sensilla, and pleopod II endopodite, to resolve a C. fukudai species complex initially detected using molecular tools. These characters showed significant size and shape variations across the geographically separated populations, of which the pleopod II endopodite displays a high taxonomic potential for yielding the new species establishment. Our pioneering attempt consequently provided the first insight into the fine-scale diversity challenging the subsequent cryptic species problem in the microcerbrid taxonomy. Coxicerberus anfidicus Messana, Argano, and Baldari, 1978 complex exhibiting the morphological inconsistency between specimens from Somalian and Maledive populations [44], for instance, would be the next case study of cryptic diversity. Finally, we expect the possible methodological application to the taxonomic impediment in other interstitial fauna as the homogenous morphology resulted from convergent evolution is a well-known phenomenon across the distant areas [104].
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-444 1/13/4/515/s1, Table S1: Specimens information examined in this study with accession numbers as well as GPS addresses; Table S2: Accession numbers of sequences used in this study's genetic analyses; Table S3: Centroid size of character employed for Figure 3; File S1: TPS files containing landmarks digitization of characters; File S2: Morpho J save files of the entire morphometric analyses.