Shell and appendages variability in two allopatric ostracod species seen through the light of molecular data

Ostracod crustaceans are among the most abundant microfossil animals. Understanding intra-and interspecific variability of their shell is of pivotal importance for the interpretation of paleontological data. In comparison to appendages, ostracod shell displays more intraspecific variability (in shape, size, and ornamentation), often as a response to environmental conditions. Shell variability has been studied with sophisticated methods, such as geometric morphometrics (GM), but the conspecificity of examined specimens and populations was never tested. In addition, there are no GM studies of appendages. We build on previously published high cytochrome c oxidase subunit I ( COI) divergence rates among populations of a brackish water species, Ishizakiella miurensis (Hanai, 1957). With landmark-based GM analyses of its shell and appendages, and additional genetic markers ( ITS, 28S, 18S ), we test if the genetic variability is structured in morphospace. This approach is the core of integrative taxonomy paradigm which has been proposed to bring the gap between traditional taxonomy and other disciplines such as evolutionary


Introduction
Ostracods are ecologically versatile crustaceans with a (mostly) well-calcified shell and an abundant and diverse fossil record.Their signature on the geological time scale can be traced back to the Ordovician (Siveter et al., 2014).Of the estimated 45,000 described species (Brandão et al., 2019) almost 80% are fossil.Long evolutionary history, abundance, and ecological versatility gave these microfossils an important role in a wide range of paleontology related studies (Boomer et al., 2003).Since a negligible number of fossilized ostracods have appendages preserved (Siveter et al., 2013(Siveter et al., , 2014)), understanding the shell morphology, ontogeny, chemical composition, variability, etc. is of pivotal importance for an accurate interpretation of their assemblages.In comparison to numerous characters provided by up to eight appendages, the shell is poor in details.Nevertheless, various features, from the shape to fine structures on the surface (internal and external), successfully contributed not only to description of so many species, but also to phylogenetic analyses of some groups (Hunt, 2007a) and comprehension of their modes of evolution (Hunt, 2013).
Ostracod shell often displays a high instraspecific variability (Aiello et al., 2016;Hunt 2007b;Wrozyna et al., 2016).One of the most notable examples is of Cyprideis torosa (Jones, 1850), where a strong variability in noding and morphology of pores is observed in correlation to salinity (Frenzel et al., 2012(Frenzel et al., , 2016)).Since paleontological studies provide no direct estimate of gene flow, only sequential changes in morphology coupled with temporal and geographic scales provide a proxy for evolutionary lineages.For example, Cronin (1985Cronin ( , 1987) ) studied a marine ostracod genus with extensive temporal and spatial scales, and found enough morphological characteristics to delineate several species, which evolved in a very short period of time after the rise of the Isthmus of Panama.This is, however, followed by a pronounced shell variability across populations of those species, which the author regarded as ecotypes, persisting during the last three million years.His methods involved multivariate analysis, but only of the linear morphometric data.Recently, there have been several studies that used 2D GM to quantify ostracod shell variability in correlation with environment conditions (see Wrozyna et al., 2018 and references therein).Although based on living species (species with appendages), none of these studies used genetics to test if there was congruence between GM and molecular data.Without genetic data it is hard to distinguish within and between-lineage morphological variation (Van Bocxlaer & Hunt, 2013).Some studies combine GM of the shell with the linear morphometrics of the appendages (Iepure et al., 2008), or use only linear morphometrics of the appendages (Wrozyna et al., 2019) to discern morphological variability biology.The results show that it is the shell shape, and not the shape of appendages, that mirrors the molecular phylogeny, and we describe a new species.Our results suggest that the landmark-based GM studies may be useful in paleontological datasets for closely related species delineation.We implement molecular clock and population statistics to discuss speciation processes and phylogeography of the two congeners in Korea and Japan.between geographically isolated populations of one species.However, there are only four published studies that attempted to couple shell variability and genetics, two using linear morphometrics (Karanovic, 2015;Macario-González et al., 2018) and two using GM (Karanovic et al., 2019;Koenders et al., 2017).Linear morphometrics has limitations in comparison to the GM.It consists of length, depth, and width measurements and thus provides little information about the shape of the studied structures, and perceived variability is strongly correlated with the size (Zelditch et al., 2004).Empirical research corroborates the advances of geometric over linear morphometrics (Mutanen & Pretorius, 2007;Fabre et al., 2014;Schmieder et al., 2015;Gabelaia et al., 2018 and references in there).
Human eye observations of appendages morphology indicate their conservative nature, even on the genus level (Karanovic, 2012;Meisch, 2000).However, there are no landmark-based GM analyses of shape variation to confirm this.Since the advent of molecular tools, the number of cryptic species has risen in all animal groups (Pfenninger & Schwenk, 2007) and ostracods are no exception (Bode et al., 2010;Brandão et al., 2010;Koenders et al., 2012;Martens et al., 2012;Schön et al., 2013).Although they conform with species concepts based on monophyly and absence of genetic intermediates, and the International Commission of Zoological Nomenclature allows for a new species to be described based on DNA sequences alone (providing that the tissue sample from the type specimen is designated, ICZN, 1999, Article 72.5.6), the newly discovered cryptic species mostly remain undescribed.The lack of easily perceived morphological differences between genetically distant species has been successfully tackled with integrative taxonomy, combining molecular results with landmark-based GM (Chesters et al., 2012;Davis et al., 2016;Karanovic et al., 2016Karanovic et al., , 2018)).Such analyses are missing in ostracods, despite the fact that they could be extrapolated to paleontological studies, where phenotypic variability can be more confidently translated into discrete species or regarded as an instraspecific variability (Butlin & Menozzi, 2000).We approach this problem by studying a living East Asian species with known high genetic diversity (Yamaguchi, 2000).
Ishizakiella miurensis (Hanai, 1957) is a brackish water species, living on muddy sand of river mouths in Japan (Hanai, 1957;Tsukagoshi, 1994;Yamaguchi, 2000) and Korea (Smith et al., 2014;Karanovic et al., 2017).The genus Ishizakiella McKenzie & Sudijono, 1981 accounts for six species, two of them known from the fossil (Yoo et al., 2012).Of the four living species, I. novaezealandica (Hartmann, 1982) is known from New Zealand, while the other three (I.miurensis, I. ryukyuensis Tsukagoshi, 1994, andI. supralitoralis (Schornikov, 1974)) have been reported from East Asia only.Yoo et al. (2012) provided a taxonomic key to species, and while there are sufficient morphological differences in their shell to distinguish fossil and recent species, major differences between living species are in the appendages, especially in the male copulatory organ.All living species have similar ecology to I. miurensis (see Tsukagoshi, 1994), with ovoviviparous mode of brood care, and no planktonic stages (see Yamaguchi, 2000Yamaguchi, , 2003)).Therefore, their distribution should be limited and species are highly susceptible to geographic isolation.
In an attempt to reconstruct phylogenetic relationships of living Ishizakiella species, Yamaguchi (2000) used a partial COI sequence from 29 Japanese and two New Zealand populations.The results shed some light on the evolutionary pattern of the genus, and revealed exceptional molecular diversity (intraspecific distances being close to 20%) and haplotypes endemicity in I. miurensis and I. ryukyuensis.The reconstructed tree suggested Downloaded from Brill.com 11/01/2023 02:11:37PM via Open Access.This is an open access article distributed under the terms of the CC-BY 4.0 License.https://creativecommons.org/licenses/by/4.0/some distinct clades within those two species, the former consisting of the Sea of Japan and Pacific clades, and the latter of two clades restricted to two groups of Nansei Islands (Amami/Okinawa and Ishigaki/Iriomote).Interestingly, I. supralitoralis had much lower intraspecific genetic distances (maximum 6%) than I. miurensis or I. ryukyuensis, although the sampling localities were more disjunct (northern part of Hokaido and southern part of Kyushu).The instraspecific distances within I. novaezealandica were about 10% between populations collected from North and South Islands.Yamaguchi (2000) did not study morphological diversity in any species, but Yamaguchi (2003) reported one population of I. ryukyuensis from Kyushu with a distinct morphology of the male copulatory appendage.The author did not describe a new species, mainly due to the inconclusive molecular phylogeny, which suggested that Kyushu and Amami/Okinawa populations form one clade, with Ishigaki/Iriomote population being its sister clade, albeit with low bootstrap support.
Such high genetic distances in mitochondrial COI suggest a presence of several, potentially cryptic, species.Crustacean lineages differ in the divergence rates, but Léfebure et al. (2006) proposed 16% distance as a threshold for species delineation in this phylum.However, even 3% was used to delineate brachiopods (Jeffery et al., 2011).In various decapod families the interspecific distances range from 6% to 20% (Da Silva et al., 2011).In ostracods, the COI divergence rates between well-defined morphological species of one genus vary between 5% and 14% (Karanovic, 2015;Martens et al., 2013), while apparently no morphological differences have been found in some cosmopolitan species where the average pairwise distances between allopatric populations can reach 20% (Bode et al., 2010).
The aim of this study is to test weather landmark-based GM can detect groups defined by sequenced genetic markers.Beside valves, we use antenna (second antenna) and hemipenis, traditionally regarded as important in species delimitation (Meisch, 2000;Karanovic, 2012).This would allow us to discern how variable is the shell in comparison to the appendages and which of the two would be more useful in species delimitations.This is the first such study done on ostracods and should be seen as an important progress towards more closely linking paleontology and neontology in ostracod related studies.

Collecting and taxonomic methods
Samples were collected from four river mouths (fig.1: localities 1, 2, 9, and 17): one in Korea (37.3241667N 129.2678167 E) and three in Japan (30.7262222 N 130.9955833 E; 40.8847222 N 141.3805556 E; and 35.3194222 N 139.4827028E) with a plankton net (mesh size 100 μm) pulled over the sediment surface, and immediately fixed in 99% ethyl alcohol.Other localities on the fig. 1 (3-8, 10-16, 18-20) are from the previous publication (Yamaguchi, 2000).Ostracods were sorted and dissected under Olympus SZX12 dissecting microscope.Glass slides with dissected appendages mounted in a drop of CMC-10 mounting media (Master Company, Inc., USA) were examined and photographed using a Leica DM 2500 compound microscope, equipped with NPlan objectives and Leica DFC 295 camera.All permanent slides used for imaging were prepared in the same manner: the entire soft body was dissected in one drop of the mounting media and cover slip was gently placed on top.Scanning Electron Micrographs (SEM) were taken with a Hitachi S-4700 scanning electron microscope at Eulji University (Seoul).Species identification was carried out with the aid of the taxonomic key provided by Yoo et   and the redescription of I. miurensis by Tsukagoshi (1994).Valves of all individuals used for the GM analysis have also been measured using Photoshop ruler tool.All examined material is deposited at the National Institute of Biological Resources (NIBR) in Seoul.

Molecular methods
DNA was extracted from all studied specimens with lysis buffer that was prepared according to Williams et al. (1992).All PCR reactions were carried out in 25 μl volume (5 μl of DNA template, 2.5 μl of 10× ExTaq Buffer, 0.25 μl of TaKaRa Ex Taq (5 units/ μl), 2 μl of dNTP Mixture (2.5 mM each), 1 μl each primer, and 13.25 μl of distilled water).PCR protocols and primer pairs used for amplification and sequencing are listed in the supplementary table S1.Presence of the DNA amplicons was verified with 1% agarose gel electrophoresis.After purification (LaboPass PCR Purification Kit, Cosmo Genetech) the PCR products were sequenced with ABI automated capillary sequencer (Macrogen, Seoul, South Korea).All obtained sequences were visualized using Finch TV version 1.4.0 (http://www.geospiza.com/Products/finchtv.shtml) to check for the quality of signal and sites with possible low resolution, and corrected by comparing forward and reverse strands.BLAST algorithm (Altschul et al., 1990) was used to check the identity of obtained sequences.All sequences were deposited in GenBank with the following accession numbers: COI MH835193-MH835291, Internal transcribed spacer (ITS) MK243379-MK243396, 28S rRNA (em region) MK240519-MK240530; 28S rRNA (vx region) MK240507-MK240518, and 18S rRNA MK240531-MK240534.Information on the number of sequences per locality and DNA marker is included in the supplementary table S2.

Phylogenetic methods
Sequences were aligned in MEGA 7 (Kumar et al., 2016) with ClustalW (Thompson et al., 1994) with default parameters.COI sequences were checked for potential stop codons with ORF finder on the NCBI website (https:// www.ncbi.nlm.nih.gov/orffinder/), using the invertebrate mitochondrial code.Beside newly obtained sequences, we also used 16 published sequences belonging to Ishizakiella miurensis available on GenBank (Yamaguchi, 2000(Yamaguchi, , 2003) ) (populations from the localities 3, 8, 10-16, and 18-20 in the fig.1).To compare the variability of the COI segment we studied with the variability of the entire COI region in ostracods, we aligned our fragment with three almost complete COI sequences from three ostracod species (accession numbers: AB114300 -Vargula hilgendorfii (G.M. Müller, 1890); KP063117 -Cypridopsis vidua (O.F. Müller, 1776); and AP014656 -Fabaeformiscandona kushiroensis Hiruta & Hiruta, 2015).The K2P model (Kimura, 1980) was used to calculate pair-wise distances between individual sequences and within and between group distances for all markers.Groups were defined based on the results of the phylogenetic analyses: sequences from the localities 1, 2, 4-7 formed one monophyletic clade and from 3, 8-20 another monophyletic clade.For the best fit evolutionary model, program jModelTest 2.1.6(Darriba et al., 2010;Guindon & Gascuel, 2003) was used with the Akaike information criterion (Hurvich & Tsai, 1989), which proposed HKY+G model of molecular evolution for the COI dataset (Hasegawa et al., 1985) and F81 for the ITS (Felsenstein, 1981).Likelihood ratio test for deviation from molecular clock implemented in DAMBE5 (Xia, 2013) accepted the molecular clock for the COI dataset.Bayesian Inference, implemented in BEAST v2.5 (Bouckaert et al., 2014) set to the 0.02, an average value for the commonly used rates of the mitochondrial COI mutations in crustaceans of 1.4% and 2.6% per million years (Knowlton et al., 1993).Yule process (Gernhard, 2008) was used as a tree prior, with BEAST default log normal distribution of the species birth rate.The analysis run for 10,000,000 generations, sampling every 1,000 generations.Software Tracer (Rambaut et al., 2018) was used to visualize results of the BEAST analyses and FigTree v1.4.3 for tree visualization.
All population genetic test statistics were calculated on the COI alignment of the four newly sampled populations (1, 2, 9, and 17) implemented in Arlequin 3.5.2.2 (Excoffier & Lischer, 2010).For each population we calculated two measures of DNA polymorphism: nucleotide and haplotype diversity indices (Nei, 1987).To detect potential geographic barriers to gene flow between populations across the studied area, F-statistics (Fst) (Hudson et al., 1992) and statistical significance of the values were computed.Demographic changes (expansion, bottleneck) in populations were estimated based on the correlation between segregating/polymorphic sites and the average number of pairwise differences (Tajima's D test) (Tajima, 1989), and using coalescent model such as Fu's test statistics (Fu, 1997).Statistical significance of Tajima's D and Fs tests was evaluated with 10,000 permutations.Mantel test (Mantel, 1967), implemented in PASSaGE 2 (Rosenberg & Anderson, 2011), was used to determine the relationship between genetic distances (K2P distances) and geographic distance to test isolation-by-distance model using 1,000 permutations.The geographic distances were calculated using great-circle distance formula.Mantel test was performed on separate datasets for the two clades recovered by the molecular phylogeny.In order to calculate distances between locality 8 and localities 3, 9-20 (all belonging to one clade) following the coast and not across the land, we added one additional point at the northern end of Honshu Island (41.498554N 140.909810E), so that all distances between the locality 8 and the rest of this clade's localities were calculated via this point.Phylogenetic relationship between detected haplotypes was reconstructed with a haplotype network using Minimum Spanning Algorithm (Bandelt et al., 1999) implemented in PopART v. 1.7 software (Leigh & Bryant, 2015: http://popart.otago.ac.nz).
To quantify historical biogeography and ancestral areas at internal nodes we applied commonly used event-based method, dispersal-vicariance (S-DIVA) (Ronquist, 1997) implemented in the program RASP 4.1 (Yu et al., 2015).Trees reconstructed with the Bayesian phylogenetic inference (implemented in BEAST v2.5) of the alignment of the COI sequences from the localities 1, 2, 9, and 17 were used for the S-DIVA analysis.

Geometric morphometric methods
We conducted GM analyses on three independent datasets: antennas, male copulatory appendages (CA), and valves.The number of studied individuals, antennas, CA, and valves per sampling locality are listed in table 1.In order to facilitate comparisons between left valve (LV) and right valve (RV), all RVs were horizontally flipped in Adobe Photoshop CS6.
For the valves we have chosen 58, for the antenna 8, and for the CA 7 homologous features, that served as landmarks (LM) in the GM analysis (Bookstein, 1991).LMs on the valves included, among others, cuticular pores (LMs 15,(23)(24)(25)(26)(27)(28)(31)(32)(33)(34)(35)(39)(40)(41)44,(52)(53)(54)(55)and 57), point of the widest curve on the posterior valve extension (LM 5), and point of the widest curve of the anterior margin (LM 58).All other LMs on the valves were crosses between ridges, either with the valve margins or internally (fig.3).The exact position of cuticular pores and ridge crossings was determined by dividing valve surface into naturally Downloaded from Brill.com 11/01/2023 02:11:37PM via Open Access.This is an open access article distributed under the terms of the CC-BY 4.0 License.https://creativecommons.org/licenses/by/4.0/occurring cell-like formations.Close observation confirmed that the number and position of these cells was almost invariable among individuals, sexes, and sides, and each cell was given a label (from A1 to Q12) according to its position.LMs on the antenna were either segment margins (LM 1,3,6 and 8), articulation between segments (LM 5), or origin of setae (LM 2, 4, and 7) (fig.4).All antenna's LMs were only on the second endopodal segment, preventing deformation that may be caused by variability in segments relation to one another as a result of slide preparation.LMs on the CA included only non-movable parts: tip of the distal process (LM1); end of a calcified part on the distal process (LM2); beginning of the distal process from the dorsal side (LM 3); end of the largest CA lobe from the dorsal side (LM 4); centre of the loop on the ejaculatory process (LM 5); anterior tip of the ejaculatory process (LM 6); and postero-ventral end of the distal process (LM 7) (fig.4).
Landmarks were digitized in TpsDig 2 (v.2.18) software (Rohlf, 2015) as the Cartesian (raw) coordinates.In order to estimate measurement error, SEMs of valves were digitized twice, while antenna and CA were photographed twice and then each image was digitized twice.Final datasets consisted of 274 digitisations for valves, 676 for antennas, and 292 for CA.
Algorithms implemented in MorphoJ (v.1.06d) (Klingenberg, 2011) were used for all statistical analyses.None of the datasets had object symmetry.In the first step, measured individuals were subjected to a generalized Procrustes superimposition (GPA) (Rohlf & Slice, 1990;Dryden & Mardia, 1998), to remove effects of non-shape variation (position, orientation, and scale).This method estimates size (represented as centroid size, CS) and shape (represented as Procrustes shape coordinates).Datasets were analysed with one-way parametric Procrustes ANOVA  S-DIVA analysis on the BI tree reconstructed from the COI alignment from populations 1, 2, 9, and 17.Pie charts at the internal nodes represent the calculated probabilities (relative frequencies) of alternative ancestral areas (reconstructions), and only those that are associated with over 0.9 posterior probabilities are colored.The top left insert represents the probability of suggested events on each of the four nodes (A, B, C, D), and the bottom left insert is the color code for the ancestral areas.‹ to quantify the relative amount of variations as a result of measurement errors (imaging and digitisation), fluctuating asymmetry, and directional asymmetry (Klingenberg & Mc-Intyre, 1998;Klingenberg et al., 2002).In the case of directional asymmetry only those individuals with digitisation of both left and right sides were considered (default function of the software).We also used ANOVA to test for the effects of sex, species, and locality on the size and shape variations.Species groups were defined based on the molecular phylogeny analysis of the COI alignment.Regression was used to test for the effects of size on shape, i.e., allometry (Klingenberg, 2016).However, further analyses were performed only on the size uncorrected datasets.We also used regression to visualize differences in CS between sexes and localities.After computing the variance covariance matrix of the Procrustes shape coordinates we used principal component analysis (PCA) for the multivariate analysis of shape variation in all datasets averaged by individual.For the shell dataset the PCA was performed on the entire dataset, and then on datasets of females and males separately.PCA served to visualize distribution of delimitating factors such as species and sexes in morphospace associated with principal components carrying the highest eigenvalues.Mean shape for each valve, sex, and species was created by subdividing the full valve dataset according to sexes and then averaging each by side and species.Shape changes were visualized with lollipop graphs (fig.5, and supplementary fig.S1).

Molecular phylogeny
All alignments were cut to the same length prior to the analysis.included sequences from all localities, only this dataset was used in further analyses.Besides having very few variable sites, the ITS dataset produced phylogenetic trees with very low support (posterior probability 0.35; tree not presented) for basal nodes, rendering unresolved relationship between sequences from newly sampled localities.Due to a very limited numbers of molecular differences in 28S and 18S fragments, the phylogenetic tree for those two markers were not reconstructed.The alignment built to compare the variability of the COI fragment used in our study with the overall variability of the almost entire COI region was 1444 bp long.The K2P distances ranged from 27% between Cypridopsis vidua and Fabaeformiscandona kushiroensis to 63% and 66% between Vargula hilgendorfi and C. vidua and F. kushiroensis respectively.In comparison to the almost entire COI region, the fragment used in the analysis of the Ishizakiella miurensis complex (from the position 98 to 352) was more conservative; distances were 23% between C. vidua and F. kushiroensis and 48% and 45% between those two respectively and V. hilgendorfi.
Tracer analysis of the BEAST results showed that the effective sample size for all measured parameters (posterior, likelihood, priors, tree likelihood, tree height, Yule model, birth rate, etc.) was far above the recommended 200, suggesting a sound estimation of the posterior distribution.The BEAST analyses produced a phylogram that divided studied localities in two clades with very high posterior probabilities (fig.1): one consisting of populations from the localities 1, 2, 4-7, and the other of the rest of populations.The latter clade also had a high support for the sequences from the localities 3 and 8.While all sequences collected from the locality 17 belonged to one clade, those from the locality 9 clustered in two clades.Strict molecular clock placed the time of the most recent common ancestor of all populations between 3 and 1.7 million years ago.The divergence time between the Korean population (locality 1) and populations from Japan belonging to the same clade was estimated to be between 1.5 and 0.7 million years ago.Similar divergence time was proposed between sequence from the locality 3 and those from localities 9-20, while older time (1.3 and 2 million years ago) was estimated for the divergence between the sequence from the locality 8 and the rest of the populations from the same clade.
A total of 38 unique haplotypes (fig. 1) were detected.Only localities 15 and 16, as well as 13 and 18, shared a haplotype.The highest haplotype diversity was detected at the locality 2 (0.9), moderate at the localities 9 and 17, and low at the locality 1 (0.45) (table 2).All populations had low nucleotide diversity.Tajima's D test and Fu's FS test were not consistent across populations, but their values were not statistically significant and cannot be used for further population dynamics discussions.On the other hand, between populations Fst values were very high and all statistically significant.Mantel test for the clade consisting of the populations 3, 8-20 returned a relatively high correlation between K2P and geographical distances (0.61), with statistically significant p-value (two-tailed p = 0.001).This correlation was much higher (0.91) and statistically significant on the same level for the clade consisting of populations 1,2, 4-7.Since from all other localities we had only one sequence each, population dynamic indices in those localities could not be estimated.
The S-DIVA analysis found 43 dispersals, 16 extinctions, and 13 vicariance events.Dispersal and extinction events were mostly restricted to terminal nodes, while vicariance events were restricted to basal nodes.However, posterior probabilities of terminal nodes were very low and therefore the probability of vicariance/dispersal/extinction events was very small.Only nodes with posterior probabilities over 0.9 were highlighted in fig. 2. It seems that vicariance was the major cause of the disjunct distribution and the primary mechanism for speciation in I. miurensis populations.Vicariance has a high probability at the root (A), the ancestral node of populations 1 and 2 (D), and at the ancestral node of populations 9 and 17.Node C supports a combination of dispersal and variance events between localities 9 and 17.

Geometric morphometrics analysis
ANOVA analyses indicated that the measurement error was negligible in all datasets for both shape and size (table 3).Fluctuating asymmetry was significantly larger than the individual variability in all datasets, for both shape and size components.Directional asymmetry was significantly (6.27 times) higher than individual variability only in the shell dataset and in the shape component of variation.
Of the additional effects in the valves dataset, sex was by far the greatest contributor to both shape and size variations, being respectively 44 and 62 times larger than the individual variability.It was followed by species (5 (shape) and 19 (size) times larger than the individual variability), and locality (almost 3 (shape) and 8 (size) times).In the antenna dataset sex contributed insignificantly to the overall variation of the shape and size.On the other hand, variation as a result of species and locality allocation significantly contributed to the overall variability but the F value was small for both (1.3) in the shape component, and considerably higher in the size component (8 for the species and 46 for the locality).In the CA dataset only locality significantly contributed to the overall variability in the size component.
Allometry did not have a significant influence on the shape variations; although the p-value was statistically significant in all three datasets, the predicted value for the valve dataset was 1.5%, while it was 7% for both antenna and CA.
PCA results were strikingly different between the valve dataset on one side (fig.3) and antenna and CA on the other (fig.4).The shape analysis of the valve dataset containing both sexes (fig.3) showed clear clustering in the morphospace described by the first two eigenvectors.The first eigenvector (carrying 43% of the total variability) separated sexes, with males being in the positive and females in the negative part.The second eigenvector (amounting to 8% of the total variability) separated the two monophyletic clades from our phylogenetic analysis.The center of gravity for the clade from the localities 9 and 17 was in the positive part, and for the clade from the localities 1 and 2 in the negative part of the PC2 axis.The shape analysis carried out on the datasets separated by sexes highlighted the differences between clades, with the first two eigenvectors describing 34% of the variability in the female dataset (fig.3), and 32% in the male dataset (supplementary fig.S1).
Regression analysis of CS showed clear contrasts between sexes and also between clades (see graph in supplementary fig.S1).Females from the localities 1 and 2 had larger CS than females from the localities 9 and 17.Males from the localities 1 and 2 were larger than males from the locality 17, but equal to slightly smaller than males from the locality 9.
The first two eigenvectors in the antenna dataset explained 57% of the total variability (fig.4).However, in the morphospace defined by these eigenvectors there was a complete overlap between clades and sexes The first two eigenvectors explained 69% of the total variability in the CA dataset, but there was no separation in morphospace by neither of the classifiers (fig.4).Both the CA and the antenna datasets had a few outliers, but closer examination of images of those individuals could not reveal any deformation, and even when we removed them from datasets the clustering in morphospace (as well as the results of ANOVA) did not change.
Paratypes: 5 females and 8 males, valves mounted on an SEM stub, appendages dissected on one slide each from the type locality and 12 females and 10 males with valves mounted on an SEM stub and appendages dissected on one slide each from: River Mouth, Tanegashima, Japan, 30°43'34"N 130°59'44"E, collector Hyunsu Yoo, 31/05/2012.Etymology: The species name is a Latin adjective "occidentalis", meaning western, and referring to the geographic distribution of the new species in relation to Ishizakiella miurensis.Differential diagnosis (fig.5 and supplementary fig.S1): Valves asymmetrical, with most prominent differences in the posterior (LM 5 displaced outwardly on the LV, inwardly on the RV) and dorsal regions (LMs 37 and 38 displaced antero-dorsally on the LV, antero-ventrally on the RV).Several groups of LMs on the valves define the major shape differences between females of the new species and I. miurensis.Posterior end of the dorsal margin (associated with LM group 1-8) is more rounded in the new species than in I. miurensis (LM group 1-8 is shifted in the antero-dorsal direction in I. occidentalis and in the posteroventral in I. miurensis).Postero-ventral region (associated with LMs 11-22) is more outwardly displaced, and the antero-ventral region (LMs 43-50) is shifted towards posterior (in I. miurensis those LMs are shifted either ventrally on the LV or anteriorly on the RV).Shape differences can also be observed in the position of sensory setae (LM 23-25 and 31), which are displaced more anteriorly in I. occidentalis and more posteriorly in I. miurensis.There are also differences in the position of some other sensory setae (LM 27 and 28), but they are not so prominent.Sensory setae and ridges in the central region of the valves are more conservative and display much less shape changes.
Average L of females was 111 μm for the LV, and 110 μm for the RV (supplementary table S4).Males of the species display shape changes in the same regions (LM groups) as females, but they are much smaller.In addition, while the magnitude of the shape differences in females is almost the same on both valves, in males the LV has more prominent shifts than the RV, especially in the anteroventral and postero-ventral regions.Interestingly, the same discrepancy can be observed in I. miurensis males but here the RV is more variable than the LV in these regions.
Average L of males was 101.29 μm for the LV and 99.57μm for the RV (supplementary table S4).
Revised diagnosis for I. miurensis is given in the supplementary note.Yamaguchi (2000Yamaguchi ( , 2003) ) recognized two clades within Japanese populations of I. miurensis: the Sea of Japan clade (corresponding to the localities 4-8) and the Pacific clade (corresponding to the localities, 3, 10-16, and 18-20) based on molecular phylogeny inferred from the same fragment as in our study.We described I. occidentalis to accommodate populations from South Korea and Japan (localities 2, 4-7).This decision was based on the type locality of I. miurensis, situated in the Kanagawa Prefecture (Hanai, 1957) within the distribution range of the Pacific clade.Human eye observations of the shell line drawings and SEM-s (Hanai, 1957;Tsukagoshi, 1994;Yamaguchi, 2000Yamaguchi, , 2003, and present study) can detect a certain amount of variability in the postero-ventral region of the female valves, however, there is no regularity.This was traditionally regarded as morphological stasis in ostracod evolution, with many species having exceptional spatial (with remarkable geographic barriers between populations) and temporal (persisting for over 5 million years) ranges (see Cronin, 1985Cronin, , 1987)).Our results imply that the application of landmark-based GM can reveal consistent morphological differences between species despite intrapopulation variability.Hopkins & Lidgard (2012) analysed extensive datasets and showed that stasis should be taken with caution, because traits have variable modes of evolution.Our data exemplify this, because, unlike shell, examined traits of CA and antenna have highly unstructured variability and only localities to some extent contribute to the shape variability of the antenna.In addition, we assume that antennas are important sensory organs in both males and females and as such exhibit developmental constrains for optimal sensory functioning.Small morphological differences between the two species can be explained by their allopatric distribution, a phenomenon that has been observed in other animal groups as well (Kiyoshi & Hikida, 2012;le Port et al., 2013;Amor et al., 2014, etc).Therefore, we can expect that species in sympatry may develop morphologically different copulatory appendages as they usually provide aid to prezygotic isolation.For example, Tsukagoshi (1988) showed that in several closely related ostracod species morphological differences in CA are enlarged when their distribution ranges overlap, a phenomenon known as character displacement.Previous studies on recent ostracod shell shape variation evidenced the presence of distinct morphotypes in allopatric populations (Wrozyna et al., 2016(Wrozyna et al., , 2018)), but since they did not study genetic markers or soft parts, our results are at the moment unmatched.We show that the valve shape distinguishes better closely related species than appendages.In addition, previous studies suggested that using structures on the shell surface alone, or in combination with several marginal points as LMs, can yield more reliable results in a GM analysis of the shape variation than outline alone (Karanovic et 2017,2018).One of the downsides of the LMs on the shell surface (such as cuticular pores) is that they may be difficult to homologize across variety of taxa.This is not the case in closely related species, because such structures tend to remain conservative and easy to locate.Another potential problem is with species with poor or no ornamentation at all, where outline based GM would be easier to apply (see Karanovic et al., 2017b).

Discussion
Although the COI sequence we analysed was relatively short, it appears that this segment is less variable than the entire COI region in ostracods, which further amplifies the divergence value between Ishizakiella miurensis and I. occidentalis.In both species there is a significant correlation between genetic and geographic distances of populations, indicating that the distribution is more likely a result of an active migration, not a passive transport.Yamaguchi (2000) suggested that I. miurensis distribution is a result of repeated colonizations of the Japanese archipelago from the Korean peninsula via land bridges during the Late Pleistocene, and the one occurring during the last glacial maximum was a distribution pathway for the Sea of Japan clade.Our S-DIVA analysis support that I. miurensis and I. occidentalis evolved by vicariance, with some dispersal events between populations of the former species.According to the very high, statistically significant fixation index (Fst value), all genetic variations between studied populations (1, 2, 9, 17) can be explained by populations structure, with all populations being highly differentiated.However, reproductive isolation could not be confirmed and therefore we base our decision to describe a new species on the phylogenetic and morphological species concepts.With statistically non-significant Tajima's and Fu's test, indicators of each populations demographic changes (bottleneck and/or expansion) could only be drawn from the haplotype and nucleotide diversities.A clear signature of a rapid growth from a small population (expansion), could be detected in the population 2 of I. occidentalis which had a high haplotype (Hd) and a low nucleotide diversities (π).On the other hand, the Korean population of I. occidentalis, with much smaller Hd and a small π, supposedly went through a bottleneck (Avise et al., 1987;Avise, 2000).For the two I. miurensis populations we included in the population statistic analysis (9 and 17), the Hd and π are suggesting expansion as well, but not as strongly as in the case of I. occidentalis from the population 2.
We agree with Yamaguchi (2000) that the present day distribution is a result of the postglacial faunal shift, but disagree on the importance of the land bridges as migration routes.Beside the fact that general trend of faunal movement during interglacial periods has a south-north trend in the northern hemisphere (Pielou, 1991), recent studies indicate that during the last glacial maximum Korea and Japan remained separated with a narrow (20 km) and shallow (10 m deep) channel (Park et al., 2000).The biogeography of the Sea of Japan clade, as defined by Yamaguchi (2000), can be explained either by spreading of the populations from some glacial refugia in the southern part of the Sea of Japan, or by subsequent northward colonization following sea transgression and reduction of the coastline.South-north trend also follows opening of the warm Tsushima current.Polyphyletic nature of the Sea of Japan clade refutes the first scenario.We speculate that the ancestors of I. miurensis inhabited southern tips of the Japanese archipelago and spread northward along the Pacific coast, following the warm Kuroshio Current.The phylogenetic tree reconstructed from the COI alignment, as well as high K2P distances, also suggest that populations 8 and 3 of I. miurensis may represent separate species, but no morphological data Downloaded from Brill.com 11/01/2023 02:11:37PM via Open Access.This is an open access article distributed under the terms of the CC-BY 4.0 License.https://creativecommons.org/licenses/by/4.0/and only one was available from each locality, preventing further systematic changes.
Frequency of splitting, level of withinspecies lineages persistence, and the length of speciation duration are the factors influencing the speciation rate, which can be interpreted from a phylogeny (Dynesius & Jansson, 2013).Our molecular phylogeny with several population splits and long branches suggests a long speciation, with high persistence of withinspecies lineages and infrequent splitting.A relatively complex genetic subdivision, with high divergence rates between some populations of I. miurensis supports this as well.All those factors are under strong influence of the lineage biology and environment.Dynesius & Jansson (2013) suggested that a low dispersal ability has a positive influence on the withinlineage persistence.Ishizakiella species do not have planktonic larvae and have a very narrow ecological niche, therefore their dispersal ability should be low, here supported by high Fst values.Of the environmental factors, spatial barriers and climatic variability should enhance splitting and speciation rates.Our molecular time tree suggests that the most recent common ancestor of the two species lived at the beginning of the Pleistocene epoch -the latest period of repeated glaciations resulting in the sea cooling and the level dropping, creating land barriers for marine organisms' dispersal.According to Cronin (1985Cronin ( , 1987) ) the frequency and duration of climatic events have stronger influence on the ostracod evolution than the magnitude of climatic changes, i.e., that frequent changes in climate, such as Milanković's cycles, result in unchanged species, or cause the morphological stasis.
Our paper has a direct implication for systematics of ostracods and beyond as formal naming of species is important for species management and possible conservation polices (Mace, 2004;Delić et al., 2017).Given that ostracods are among the most abundant organisms in the fossil record, this study emphasizes the importance of explicit quantification of shape through GM which opens a new perspective to paleontological studies.

Figure 1
Figure 1 Biogeography and molecular phylogeny of Ishizakiella miurensis (Hanai, 1957) and I. occidentalis sp.nov.The map shows newly sampled localities (numbers in squares) and localities from Yamaguchi (2000), while the line across Kyushu and Honshu represents proposed distribution limits of the two species.Locality numbers correspond to the numbers on the haplotype network (insert between Korea and Japan) and on the Bayesian inference tree constructed from the COI dataset (at the bottom).Numbers above branches on the tree represent posterior probabilities; only values above 0.95 are presented.Numbers between brackets are 95% confidence intervals of the divergence times (in millions of years) of the main clades, based on the strict molecular clock.The scale bar shows the estimated number of substitutions per site.Downloaded from Brill.com 11/01/2023 02:11:37PM via Open Access.This is an open access article distributed under the terms of the CC-BY 4.0 License.https://creativecommons.org/licenses/by/4.0/ Figure2S-DIVA analysis on the BI tree reconstructed from the COI alignment from populations 1, 2, 9, and 17.Pie charts at the internal nodes represent the calculated probabilities (relative frequencies) of alternative ancestral areas (reconstructions), and only those that are associated with over 0.9 posterior probabilities are colored.The top left insert represents the probability of suggested events on each of the four nodes (A, B, C, D), and the bottom left insert is the color code for the ancestral areas.
Figure 3 Landmark-based GM shape analyses of the valves.Target symbols with numbers on the LV SEM photograph of Ishizakiella occidentalis sp.nov., female, holotype (top) represent landmarks.The left scatter plot shows delimitation of species and sexes, and the right scatter plot shows delimitation of females only, based on the landmarks for specimens collected from localities 1, 2, 9, and 17.Ellipses represent 95% confidence intervals, and numbers between brackets represent percentages of total variance for the first two eigenvectors (PCs).

Figure 4
Figure 4 Landmark-based GM shape analyses of the antenna (upper left, Ishizakiella occidentalis sp.nov., female, holotype) and CA (upper right, I. occidentalia sp.nov., male, allotype).Target symbols with numbers on the light photographs of the appendages represent landmarks.Scatter plots below the appendages are respective graphical visualizations of the principal component analysis based on the landmarks for specimens collected from localities 1, 2, 9, and 17; ellipses represent 95% confidence intervals, and numbers between brackets represent percentages of total variance for the first two eigenvectors (PCs).

Figure 5
Figure 5 Female shape changes (scaled 5×) associated with species and valves relative to the mean shape of the sample consisting of females only.Landmark numbers as in the fig.3.

Table 1
Sampling localities with numbers of individuals tested for genetic markers and morphological variations

Table 2
Pairwise Fst values and population statistics based on the COI alignment of four populations(1, 2, 9,  and 17) Downloaded from Brill.com 11/01/2023 02:11:37PM via Open Access.This is an open access article distributed under the terms of the CC-BY 4.0 License.https://creativecommons.org/licenses/by/4.0/

Table 3
Shape and size variation inferred by a one-way parametric ANOVA.Only Goodall's F critical values (F) and probability of finding a random value larger than the observed (P) are shown.Imaging error was only measured for antenna and CA Downloaded from Brill.com 11/01/2023 02:11:37PM via Open Access.This is an open access article distributed under the terms of the CC-BY 4.0 License.https://creativecommons.org/licenses/by/4.0/