Geological and paleoclimatic events reflected in phylogeographic patterns of intertidal arthropods (Acari, Oribatida, Selenoribatidae) from southern Japanese islands

A comprehensive study of the intertidal oribatid mite fauna of southern Japanese islands revealed the presence of the selenoribatid Arotrobates granulatus Luxton, 1992 and two yet undescribed species. The latter are herein described as Indopacifica taiyo n. sp., occurring from the Southern to the Central Ryukyus, and Indopacifica tyida n. sp., which was only found on the most western island of the Ryukyus, namely Yonaguni. A concomitant molecular genetic study using mitochondrial COI and 18S rRNA gene sequences, demonstrated that the phylogeographic pattern of I. taiyo n. sp. reflects recent expansion on the Southern and Central Ryukyus, probably due to existing land bridges during the late Pleistocene. Arotrobates granulatus , on the other hand, shows three distinct lineages, one on Japanese mainland, another on the island of Amami, and the third on part of the Central and Southern Ryukyus. These lineages are most likely the result of the break- up of a large peninsula reaching from China to the Northern Ryukyus about 1.2– 1.7 million years ago. Despite emerging land bridges in the late Pleistocene, this species was not able to expand its range again which indicates very low dispersal abilities. Morphometric data of I. taiyo n. sp. show considerable intraspecific variation between island populations correlating with geography. This found variation is suggested to be a result of phenotypic plasticity caused by diverging local environmental factors. From an ecological perspective, all three found species are classified as intertidal rock- dwellers feeding on diverse algae, whereas I. taiyo n. sp. and Arotrobates granulatus occasionally occur in mangrove habitats.

The Selenoribatidae are the most diverse group among these marine-associated oribatid mites, with nine genera and 32 species worldwide. Most members of this family were reported from the Indo-Pacific area (Pfingstl & Schuster, 2014;Procheş & Marshall, | 3 PFINGSTL eT aL. taxa, and the fauna is largely divided among the above-mentioned three geographic areas which is supposed to be partly caused by these geographic gaps (Muraji et al., 2012).
A recent phylogeographic study on intertidal mites from the family Fortuyniidae (Pfingstl, Wagner, et al., 2019) showed a slightly divergent picture. There were no endemic species and while the Tokara gap was strongly reflected in the morphological and molecular genetic sequence data, indications of the Kerama gap acting as an effective barrier to dispersal were completely lacking. This pattern was suggested to be the result of low sea level during the last Pleistocene (Ujiié et al., 2003) which allowed recent expansion and gene flow between several island populations of fortuyniid mites (Pfingstl, Wagner, et al., 2019). Apparently, these tiny wingless organisms show better dispersal abilities than partially volant terrestrial animals, such as insects and birds, and it was suggested that their ability to drift on water may be a key factor (Pfingstl, Wagner, et al., 2019).
During the field work for this latter study, we were able to additionally collect numerous populations of selenoribatid mites from various Japanese Islands, which gave us the opportunity to compare their phylogeographic patterns with those of the abovementioned closely related Fortuyniidae. Therefore, aims of the present study were (a) to determine collected taxa and describe possible new species, (b) to infer phylogeographic scenarios for each found selenoribatid group and compare these to patterns of fortuyniid mites, (c) to interpret the found similarities or divergences, and (d) to update biogeographic information for Japanese intertidal mites.

| Sample locations
Samples of littoral algae were scraped off rocks, concrete walls, and mangrove roots with a small shovel and then put in Berlese-Tullgren funnels for 12-24 h to extract mites. Afterward, specimens were fixed in ethanol (100%) for morphological and molecular genetic investigation. The majority of samples was taken by Shimano, Hiruta, and Pfingstl and if not, the name of the collector is provided.
Each sample location was given a code; these codes are presented in parentheses and will be used throughout the manuscript to allow easier linking of geographic information. The suffixes "-jima/-shima" and "-gawa" are Romanized Japanese meaning "island" and "river," respectively; in the graphs and tables showing results, these suffixes are not given due to a shortage of space.

| Molecular genetic analyses
For this study, 121 specimens of the family Selenoribatidae (67 of genus Indopacifica, 53 of the genus Arotrobates and one outgroup specimen of T. barbara) were analyzed. Whole genomic DNA was extracted from ethanol-fixed individuals using Chelex resin according to the adjusted protocols in Pfingstl, Lienhard, et al. (2019). We amplified a partial sequence of the mitochondrial DNA cytochrome c oxidase subunit I (COI) gene (600 bp amplicon length) and the complete 18S ribosomal RNA gene (ca. 1900 bp amplicon length) of major mitochondrial lineages, following the protocols of Pfingstl, Lienhard, et al. (2019). PCR conditions and primers as well as sources are summarized for each locus in Table S1. Subsequent DNA purification steps included enzymatic ExoSAPIT (Affymetrix) and Sephadex G-50 resin (GE Healthcare). Cycle sequencing, using BigDye Sequence Terminator v3.1 kit (Applied Biosystems), was conducted according to Schäffer et al. (2008). Automatic capillary sequencing and sequence visualization was operated on an ABI3500XL (Applied Biosystems) device. Furthermore, already published sequences of I. parva (18S rRNA: MH285690 and COI: MH285671) and I. pantai (18S rRNA: MH285691, MH285692 and COI: MH285654, MH285652) (Pfingstl, Lienhard, et al., 2019) were downloaded from GenBank and included in the alignment.
MUSCLE (Edgar, 2004), integrated in the software MEGA 7.0 (Kumar et al., 2016), was used to align sequences. Concatenation of the single locus datasets and indel coding of the 18S alignment was performed in the perl-based program 2matrix (Salinas & Little, 2014). Finally, three alignments were created and analyzed (a) COI (564 bp), (b) 18S (1964 bp; incl. indels and binary coded gaps), and (c) COI + 18S (2528 bp; incl. indels and binary coded gaps). For more detailed information, see alignments in Data S1S-3. All generated sequences have been deposited at GenBank and are accessible at following numbers, MW289085-MW289203 (COI) and The best fitting model of molecular evolution was determined based on Bayesian Information Criterion (BIC) in PartitionFinder v2.1.1 (Lanfear et al., 2017). Additionally, for the COI dataset, based on all three codon positions as starting scheme, a "greedy" partition search (Lanfear et al., 2012) in PartitionFinder v2.1.1 was conducted. The results of these analyses are summarized in Table   S1. Phylogenetic inference was conducted for each single locus and the concatenated dataset in RaxML HPC v.8.0. (Stamatakis, 2014) and MrBayes v3.2 (Ronquist et al., 2012). Node support for the maximum likelihood framework was assessed using 10.000 bootstrap replicates and GTR-GAMMA model of evolution.
Posterior probabilities for the Bayesian inference were obtained from Metropolis-coupled Markov chain Monte Carlo (MCMC) simulations using two independent runs, eight chains, a sampling frequency of 10,000 and the respective model of evolution. The analyses were run for 20e −6 , 15e −6 , and 10e −6 generations for the concatenated, COI and 18S dataset, respectively. To ensure parameter convergence and stationarity of chains, split deviation frequencies lower than 0.01 were assured and effective sample sizes (ESS) were investigated in Tracer 1.6  (available at http://beast.bio.ed.ac.uk/Tracer). The first 25% of obtained trees of the MrBayes run were discarded and the rest summarized in terms of a 50% majority rule consensus tree. All trees were visualized in FigTree v1.4.2 (Rambaut, 2014) (available at http://tree.bio.ed.ac.uk/softw are/figtree). Net interspecific pdistances between groups, using 1000 bootstraps replications, were calculated in MEGA 7.0. Minimum interspecific and the maximum intraspecific divergence (in %) were calculated in R vs. 3.6.1 (R Core Team, 2017) using the functions "nonConDist" and "max-InDist" from the R-package SPIDER v.1.5 (Brown et al., 2012). Raw intraspecific genetic distances were obtained for both genera by the function "dist.dna" implemented in the R-package APE v.5.4 (Paradis et al., 2004) and visualized using the package PHEATMAP v.1.0.8 (Kolde & Kolde, 2015). Furthermore, to infer the phylogeographic relationships among genera, TCS networks (Clement et al., 2002) were calculated in PopART (Leigh & Bryant, 2015; available at http://popart.otago.ac.nz).

| Morphometric analyses
Specimens were embedded in lactic acid for temporary slides, and measurements were made using a compound light microscope (Olympus BH-2) and ocular micrometer. A set of 15 continuous variables (bl-body length, dPtI-distance between pedotecta I, db i -distance between inner borders of bothridia, db o -distance between outer borders of bothridia, nw da -notogastral width on level of seta da, nw dm -notogastral width on level of seta dm, nw dp -notogastral width on level of seta dp; cl-camerostome length, cw-camerostome width, dcg-distance between camerostome and genital orifice, dac3-distance between acetabula 3, glgenital orifice length, gw-genital orifice width, al-anal opening length, aw-anal opening width, see Figure S1) were measured in 75 Indopacifica specimens from six different populations originating from four different Ryukyu islands (Yonaguni-jima, Iriomotejima, Ishigaki-jima, and Okinawa-jima). For species discrimination, all populations from a species were pooled; for the comparison of populations in the wider distributed Indopacifica taiyo n. sp., 61 specimens from five different locations from above given four islands were used for analysis. Specimens used for morphometric comparison were not the same as used for molecular genetic analyses, but they belonged to the exact same populations (same sample, approx. 10 cm 2 patch of alga). Specimens were sexed based on internal structures (spermatopositor, ovipositor), and sexes were separated for morphometric analyses.
For univariate statistics of Indopacifica, minimum, maximum, mean, and standard deviation for each variable were calculated.
Multivariate analyses investigating differences between putative Indopacifica species included a Principal Component Analysis (PCA; using a variance-covariance matrix) and a Non-metric Multidimensional Scaling (NMDS; based on Euclidian distances, twodimensional); both analyses were performed on log 10 -transformed size-corrected data. No rotation was applied to the multivariate data. Size correction was done by dividing each variable through the geometric mean of the respective specimen.
For the investigation of divergence among different populations of Indopacifica, a Kruskal-Wallis and a Mann-Whitney U-test (Bonferroni-corrected p-values) were used for comparing the means of variables for pairwise comparisons to clarify whether single variables differ significantly between the populations. Additionally, a Discriminant Analysis (LDA) was performed on log 10 -transformed raw data to show size and shape differences between populations.
All analyses were performed with PAST version 3.11 (Hammer et al., 2001).
Due to the relatively low number of A. granulatus specimens available, we were not able to perform a morphometric analysis of this species.

| Drawings and photographs
Preserved specimens were embedded in Berlese mountant for microscopic analysis in transmitted light. Drawings were performed with an Olympus BH-2 Microscope equipped with a drawing attachment, and then they were scanned and afterwards processed and digitized with the free and open-source vector graphics editor Inkscape (https://inksc ape.org).
For photographic documentation, specimens were air-dried and photographed with a Keyence VHX-5000 digital microscope.
Morphological terminology used in this paper follows that of Grandjean (1953) and Norton and Behan-Pelletier (2009). Etymology. The sound of specific epithet "taiyo" has the meaning of both "big ocean" and "the sun" in the general Japanese and "big ocean" refers to the wider coastal distribution from the Southern to the Central Ryukyus. It is given as noun in apposition. Gnathosoma. Chelicera chelate, with two teeth on each digit.
Distal part of rutellum developed as thin, triangular, slightly inwardly curved membrane with longitudinal incision. Setae a and m long, smooth. Mentum regular, finely granular, seta h simple, long.
Distribution. The distribution of this species ranges from the Southern Ryukyus (Yonaguni-jima, Iriomote-jima, Ishigaki-jima) to the Central Ryukyus, whereas Okinawa-jima represents its northern limit ( Figure 3). Gnathosoma. Chelicera chelate, with two teeth on each digit.
Distal part of rutellum developed as thin, triangular, slightly inwardly curved membrane with longitudinal incision. Setae a and m long, smooth. Mentum regular, finely granular, seta h simple, long.
Distribution. This species was only found in a single location on Yonaguni the south westernmost island of Japan ( Figure 3).
Remarks. Indopacifica tyida n. sp. can be distinguished from I. taiyo

| Remarks on the morphology of Arotrobates granulatus Luxton, 1992
The populations of A. granulatus from Japanese mainland ( Figure 5) and the Central and Southern Ryukyus show conformity in their morphology and no distinct differences could be detected (Table 2). A slight deviation is shown only in the specimens from Kagoshima (JP_86) which possess very fine adanal setae instead of all being reduced to their alveoli. The specimens from Hong Kong described by Luxton (1992) differ from the herein investigated Japanese individuals by lacking lamellar setae (le) and one pair of notogastral seta and by having only alveolar interlamellar (in) and epimeral setae (3a, 4a). The specimens TA B L E 1 Leg setation in detail, identical in both Japanese Indopacifica species  Karasawa and Aoki (2005) also lack the lamellar setae (le) and show only alveolar interlamellar setae (in), but they lack two pairs of notogastral setae and show variable adanal setation (Table 2). However, most of these differences are subtle and may be results of either intraspecific variation or inaccurate observation. In the case of adanal setae, Karasawa and Aoki (2005) already mentioned that they can be present or absent and this variation was also found in the present study which confirms intraspecific variation at least in some aspects. The lamellar setae which were not detected by the other authors (Karasawa & Aoki, 2005;Luxton, 1992) are very delicate structures and partly covered by cerotegument in most of the presently investigated populations and we also had initial difficulties to find them. Consequently, it is possible that these setae were present in the other studies as well but were overlooked due to their very fine nature. Moreover, the interlamellar setae in our specimens are very short and may also be easily mistaken for simple alveoli.
The differences in notogastral setation, on the other hand, may represent real divergences though there are certain discrepancies. Seta c 3 of the herein studied specimens is located in a very lateral position so that it is difficult to detect from a dorsal view. However, Karasawa and Aoki (2005)  To conclude, the diagnostic value of the existing morphological differences between specimens from the different studies is not clear and needs further investigation.

| Genetic data
The topology of the Bayesian Inference (BI) tree based on con-  Note: Japanese mainland: Omishima, Kagoshima; Central Ryukyus: Amami, Okinawa; Southern Ryukyus: Ishigaki and Iriomote. dto = difficult to observe; ? = information not given clearly in literature; c 3 /la = between seta c 3 and la; la/lm = between setae la and lm. and the Thai I. parva is placed as the sister group to I. taiyo n. sp. The latter consists of two distinct groupings, one including specimens from Iriomote-jima (JP_46, JP_50) and Yonaguni-jima (JP_56) and the other includes only specimens from Iriomote-jima (JP_47) ( Figure 6).
Indopacifica tyida n. sp., on the other hand, is represented by a homogeneous single clade. Another distinct clade is represented by specimens from Shikoku (JP_20, Japanese mainland) whereas these are placed clearly outside the Indopacifica group and thus they are given as unidentified selenoribatid genus (no specimens or useable remains of these were available for morphological identification).
Maximum intraspecific distances in the COI gene sequence are 0% in I. tyida n. sp. (single haplotype), 8.5% in I. taiyo n. sp., and 14.8% in A. granulatus (Figure 7). Mean net distances between species range from 9.8% to 14.4%, with a net divergence of 12.7% between the two Japanese Indopacifica species (Table 3).
TCS analysis resulted in a network of 13 haplotypes for I. taiyo n. sp. and only one haplotype for I. tyida n. sp. The I. taiyo n. sp.
LDA on raw data of only male I. taiyo n. sp. shows the same pattern as with all individuals but with clearer separation of each island

| Taxonomy and systematics
The genus Indopacifica was only recently erected (Pfingstl, Lienhard, et al., 2019), and with the discovery of the two Japanese species, it already comprises six different species. The morphology of this genus is quite homogeneous with very few characters allowing to distinguish between species. A similar phenomenon was reported in the intertidal genus Fortuynia where it was suggested to be a result of limited habitat preferences (Pfingstl, 2015)  Arotrobates granulatus was originally found at the coast of Hong Kong (Luxton, 1992) and later it was reported from the Southern and Central Ryukyus (Karasawa & Aoki, 2005). In this study, we could confirm the distribution on the Ryukyus and we could demonstrate additional occurrences on the Japanese mainland. As outlined earlier, there are slight morphological differences between the specimens from each study whereas most of these very few divergences may actually be the result of inaccurate observations or intraspecific variation. Accordingly, it is not possible to tell with confidence if specimens from the different studies belong to the exact same species. However, it would be highly unlikely that the specimens taken from the same Ryukyuan islands by Karasawa and Aoki (2005) and Considering the phylogeny of studied selenoribatid taxa, A. granulatus is placed on a long branch in the 18S trees, distant from the other groups, the outgroup included. There are two possible explanations for this excessively long branch: (a) the 18S gene evolves at a distinctly higher rate in A. granulatus, albeit this is rather unlikely because the rate must be three times higher than in the other groups; (b) the Caribbean T. barbara, which was chosen as outgroup due to its far distant geographic occurrence, is closer related to Indopacifica than Arotrobates is. This is visible in Figures S5-S7  Note: Net average uncorrected p-distance between species are shown (below diagonal). Standard error estimate(s) are given in grey above the diagonal.
F I G U R E 8 TCS Haplotype networks based on COI sequences for the three selenoribatid species. Each circle corresponds to one haplotype, its size is proportional to its frequency, and the number of mutations is indicated as hatch marks. Small black circles represent intermediate haplotypes not present in the dataset. Codes refer to sample locations as given in the material and method section; different islands are given in different colors not been investigated yet, neither on a morphological base nor on a molecular genetic base and we are the first to put this taxon into a phylogenetic context. Nevertheless, for a comprehensive phylogenetic study of Selenoribatidae more taxa from all geographic regions are needed and therefore phylogenetic relationships within this fam- ily remain yet uncertain.

| Phylogeography
The high genetic diversity of Japanese A. granulatus coincides with the geological history of this region. The Japanese mainland populations are clearly separated from all Ryukyuan populations which corresponds to the Tokara gap. This gap, which was presumably formed during the late Pliocene, is located on the Tokara tectonic strait and is suggested to have served since then as a barrier for the dispersal of amphibians and reptiles (Ota, 1998). An earlier study on oribatid mites from Japan (Pfingstl, Wagner, et al., 2019) also demonstrated that the Tokara gap is clearly reflected in the genetic structure of the intertidal Fortuynia shibai Aoki, 1974, and the present data confirm this strait as an effective dispersal barrier for small intertidal mites.
The haplotype network of A. granulatus further shows a distinct Amami-ohshima clade which is clearly separated from the other Ryukyuan island populations, including Okinawa-jima.
The latter island is, together with Amami-ohshima, part of the Central Ryukyus which is a chain of several neighboring islands separated by smaller stretches of open ocean. In contrast to A. granulatus, the intertidal F. shibai as well as Fortuynia churaumi Pfingstl, Shimano & Hiruta, 2019 appear to be closely related and even share haplotypes between both islands (Pfingstl, Wagner, et al., 2019), which was suggested to be a result of emerging land bridges during the late Pleistocene (Lin et al., 2002). Apparently, A. granulatus has not been able to cross these alleged land bridges probably due to a very low dispersal ability. Furthermore, the results of the haplotype network indicate hardly any closely related haplotypes on the same island or even at the same locality. This implies extremely low gene flow even between geographically close populations and could explain why A. granulatus was not able to cross the oceanic stretch between Amami-ohsima and Okinawa-jima. The same barrier for gene flow was found in certain insects, for example the wood-feeding cockroach Salganea taiwanensis ryukyuanus Asahina, 1988(Maekawa et al., 1999 or the firefly Curtos costipennis (Gorham, 1880) (Muraji et al., 2012).
The channel between the Amami and the Okinawa region opened sequentially around 1.3-1 million years ago resulting in isolation between the islands, and despite emerging super islands and land bridges in the late Pleistocene (Muraji et al., 2012), it apparently remained an insuperable barrier for the above-mentioned taxa. An Atlantic Fortuynia species was observed to show a specific floating behavior when washed into open water which could facilitate dispersal via ocean currents (Pfingstl, 2013). Given that the Japanese Fortuynia species show the same behavior, this could explain why  163-197 (179 ± 8.3) 154-194 (174 ± 11.1) nw dm 175-203 (190 ± 5.6) 154-194 (176 ± 13) nw dp 139-182 (157 ± 10) 151-172 ( (Yamamoto et al., 2006) suggested that this peninsula began to break 1.6 million years ago. Based on this scenario and the given distribution from Hong Kong (Luxton, 1992) to Japanese mainland, it is possible that an ancestral A. granulatus population was distributed on this peninsula and that the break-up of this landmass induced a sequential vicariant process reflecting the found pat-  (Lin et al., 2002). In contrast to A. granulatus, I. taiyo n. sp. was apparently able to use these late Pleistocene low sea level stands to disperse across the Southern and Central Ryukyus. This assumption is further corroborated by the fact that two other Ryukyuan intertidal mites, F. churaumi and Alismobates reticulatus Luxton, 1992, show comparable genetic patterns which were also suggested to be a result of the same evolutionary scenario (Pfingstl, Wagner, et al., 2019).
Apart from all the closely related island populations, there a two genetically distinct haplotypes in I. taiyo n. sp., one from Okinawajima and one from Ishigaki-jima. A similar distinct clade from Okinawa-jima was found in F. churaumi, whereas it was proposed that it either represents a relic of more ancient populations or the result of a recent colonization event from more isolated areas of this region (Pfingstl, Wagner, et al., 2019). As two independent colonization events from far isolated areas are less likely, we think that the distinct haplotypes herein rather represent remnants of ancient populations that have subsisted since the late Pleistocene expansion, but certainly more data are necessary to reliably answer this question.
The Tokara gap as well as a barrier between the Amami and Okinawa region are clearly reflected in the sequence data of A. granulatus. However, the Kerama gap, located between the Central and Southern Ryukyus, is not detected in the phylogeographic pattern of any herein investigated species, despite its pivotal role as distribution barrier for insects (Tojo et al., 2017), amphibians, reptiles (Ota, 2000), birds (Tokuda, 1969), and mammals (Motokawa, 2000 that the Kerama gap is apparently not an effective zoogeographical barrier for intertidal mites and the present data clearly support this hypothesis. The current distribution ranges indicate that A. granulatus showed a wide ancestral distribution along the Asian coastline before the opening of the Tokara strait in late Pliocene (Ota, 1998)

| Morphological variation between populations
Morphometric data demonstrate that there is a considerable in-  (Pfingstl & Baumann, 2017). However, the reason for this variation remains unknown.

| Ecology
Arotrobates granulatus was first found in oyster shells on a rocky shore (Luxton, 1992), later it was also reported from bark of the mangrove Bruguiera gymnorrhiza and from algae on a mangrove forest floor (Karasawa & Aoki, 2005). In the present study, A. granulatus was found in different kinds of algae and once on barnacles and mussels which indicates a wider ecological range for this species. Based on results from the present study, A. granulatus seems to feed on a variety of algae but prefers rocky substrates (though it occasionally may occur in mangrove habitats).
Species of the genus Indopacifica were recently discovered for the first time (Pfingstl, Lienhard, et al., 2019) and therefore knowledge about their ecology is largely incomplete. One of the few existing studies stated that they seem to inhabit and feed on diverse kinds of algae growing on intertidal rocks (Resch et al., 2019) and the present study basically confirms this statement. Indopacifica tyida n.
sp. was only found in a single location which was a typical rocky shore and I. taiyo n. sp. occurred mainly in algae growing on rocky substrate, and only in three out of 21 samples they were found in mangrove associated habitats with very low specimen numbers.
Thus, similar to A. granulatus, both species apparently prefer rocky intertidal environments but may occasionally be found in other habitats too.
In eight locations I. taiyo n. sp. and A. granulatus were found syntopically but based on the present findings we are not able to tell whether they use different niches within the microhabitat, for example, feeding on different parts of the algae or showing a vertical zonation. Both Indopacifica species, on the other hand, never occurred syntopically and microscopic investigations revealed a blue shimmering gut content in nearly all I. tyida n. sp. specimens, indicating a different diet of this species. The found blue shimmering gut content may indicate that I. tyida n. sp. may use cellular blue-green algae (cyanobacteria) as food source while I. taiyo n. sp. may feed on large multicellular algae and this could explain why they were never found together. However, further gut content analyses and food preference tests are necessary to verify such a statement. Murakami; FY2020-FY2022).

DATA AVA I L A B I L I T Y S TAT E M E N T
Gene sequence data are available at GenBank-NCBI; https://www.