Contrasting phylogeographic patterns of intertidal mites (Acari, Oribatida) along the South African shoreline

The South African coast is known to harbor four different species of intertidal oribatid mites and their distribution strongly correlates with marine ecoregions. Relatively little is known about the dispersal of these organisms and how populations of different locations are connected. To test dispersal abilities and connectivity of these South African species, we performed a morphometric and molecular genetic study. COI gene sequences of two of the widely distributed South African intertidal oribatid mite species revealed clearly contrasting patterns. Halozetes capensis, which occurs in the Agulhas Ecoregion, shows distinct genetic structuring, whereas Fortuynia elamellata micromorpha, which is distributed in the Natal Ecoregion, exhibits gene flow between all populations. The paleoenvironmental history and specific ocean current pattern are suggested to be responsible for these patterns. During the last glacial maximum, the colder climate and the weakening of the Agulhas Current possibly resulted in a bottleneck in the warm-adapted F. e. micromorpha populations, but the subsequent global warming allowed the populations to expand again. The cold-adapted H. capensis populations, on the other hand, experienced no dramatic changes during this period and thus could persist in the Agulhas Ecoregion. Considering transport on ocean currents, the Agulhas Current could be further responsible for the connectivity between the Fortuynia populations. But the deflection of this current in the Agulhas Ecoregion could support the isolation of Halozetes populations. The concomitant morphometric study demonstrated morphological homogeneity among populations of Fortuynia and thus confirms strong connectivity. The Halozetes populations, on the other hand, form two different morphological groups not reflecting geography.


Introduction
South Africa's shoreline stretches over a distance of about 3650 km (Lombard, 2004) and consists of sandy beaches and rocky shores, made of sandstone, mudstone, granites, and shales (Mucina et al., 2006). The coastline is almost linear in outline and strongly wave exposed, especially in the southwest (Shillington & Harris, 1978), and the tidal regime ranges from 2 to 2.5 m during springtide to about 1 m during neap tide (Field & Griffiths, 1991). Oceanographic features of South Africa's coast are dominated by two major current systems ( Fig. 1), the cold Benguela Current flowing northwards along the Atlantic coast and the warm Agulhas Current moving southeastwards along the Indian Ocean coast (e.g., Griffiths et al., 2010). Six marine ecoregions are recognized, Southern Benguela, Agulhas, Natal, Delagoa, Southeast Atlantic, and Southwest Indian, whereas the first four include the coast, continental shelves, and shelf edge, and the latter two only comprise the bathyal zone and the abyss (e.g., Sink et al., 2012). Despite some controversies regarding the names, boundaries, and overlap zones, these biogeographic provinces have been established for several groups of seaweeds, coastal fish, and intertidal invertebrates (e.g., Emanuel et al., 1992;Field & Griffiths, 1991;Turpie et al., 2000).
Intertidal oribatid mites belong to the latter group and a recent study investigating the South African fauna of these animals (Pfingstl et al., 2021a) revealed the presence of four 1 3 species from three families and demonstrated that their distributions also coincide with marine biogeographic regions. Halozetes capensis Coetzee and Marshall (2003) (family Podacaridae) is restricted to the warm-temperate Agulhas Ecoregion, Fortuynia elamellata micromorpha Marshall and Pugh (2002) (family Fortuyniidae) and Schusteria ugraseni Marshall and Pugh (2000) (family Selenoribatidae) are confined to the warmer subtropical Natal Ecoregion, and Selenoribates divergens Pfingstl (2015) (family Selenoribatidae) was recorded from a single location in the tropical Delagoa Ecoregion. These intertidal mites represent a small group of tiny organisms that have colonized the marine littoral environment, but they have not yet completely crossed the ecological barrier between marine and terrestrial environments (Procheş & Marshall, 2001). Consequently, they are still air-breathing animals using plastron respiration to tolerate tidal inundation (e.g., Pfingstl & Krisper, 2014;Pugh et al., 1990) and they forage walking on eight legs but are not able to actively swim. These mites mainly feed on intertidal algae which they also use as substrate and many species are typical elements of rocky shore biota (e.g., Pfingstl, 2017). Specific studies on foraging behavior, feeding, and microhabitat preferences of known South African intertidal oribatid mites are lacking so far. Moreover, relatively little is known about the regional persistence of rocky shore assemblages over time scales greater than ca. 30 years and spatial scales larger than a few meters (e.g., Paine & Trimble, 2004) and this also includes intertidal oribatid mites. Due to their minute size (approx. 0.3-0.6 mm), these organisms are difficult to observe in the field and virtually nothing is known about how they migrate between patches of algae of a single location or how they disperse along the coastline over long distances. Three of the above-mentioned species, Halozetes capensis, Fortuynia e. micromorpha, and Schusteria ugraseni, are known to show distributions ranging over several hundred kilometers of coastline (Pfingstl et al., 2021a;Procheş & Marshall, 2002), and the almost linear outline of the South African shore with its distinct pattern of ocean currents offers the opportunity to provide first insights into connectivity and level of genetic structuring between populations of different locations.
For this purpose, we studied numerous populations of these three species from the Western Cape to KwaZulu-Natal by means of molecular genetics and morphometrics. Apart from connectivity, we also aimed to test for morphological divergence between far distant populations and to discuss possible ways of short-and long-distance dispersal in view of the present results.

Collection
Samples were collected during two fieldtrips in February and March 2019 and 2020 and an additional trip in October 2020. Patches of intertidal algae, barnacles, mussels, and tubeworm colonies were removed with a small shovel or a knife during low tide. Mobile, self-made, Berlese-Tullgren funnels were used in temporary field laboratories to extract the mites from the samples. Mites were then collected alive in small plastic vessels lined with plaster of Paris, cleaned from debris with a fine brush, and preserved in absolute pure ethanol. For microscopic investigation (Olympus BH-2 microscope) and morphometric measurements, specimens were embedded in temporary slides using lactic acid.
Graphs and maps were created using the Concepts software (https:// conce pts. app) and the free open-source vector graphics editor Inkscape (https:// inksc ape. org), and final processing was done with Adobe Photoshop 7.0.

Sample locations
For the sample locations ( Fig. 1), see Pfingstl et al. (2021a) for details on locality and habitat. On nearly each locality, several samples of algae (20-100 m apart) were taken, and for clear identification, sample IDs were given and these can be found below as simple codes (e.g., ZA_22). (

DNA extraction and sequencing
Whole genomic DNA was extracted from 78 ethanol-fixed individuals of Halozetes capensis and 121 ethanol-fixed individuals of Fortuynia e. micromorpha (individuals were extracted separately, they were not pooled). Therefore, Chelex resin was used, according to the adjusted protocols in Pfingstl et al. (2019). A ca. 560 bp long fragment of the cytochrome oxidase subunit I (COI) was amplified following the protocols of Pfingstl et al. (2019). PCR conditions and primers are shown in Table 1. DNA purification steps afterwards 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 were run on an ABI3500XL (Applied Biosystems) device. The sequences were then loaded into the software MEGA 7.0 (Kumar et al., 2016) and aligned manually. Due to unknown reasons, we could not sequence the COI gene  Table S1).

Population differentiations
The final alignments included 121 (579 bp) and 78 (555 bp) individual sequences for the genera Fortuynia and Halozetes, respectively. We calculated for the Fortuynia dataset the amount of genetic differentiation between locations as haplotype frequencies (FST) and mean pairwise differences (ΦST) in Arlequin v3.5.2.2 (Excoffier et al., 2005), whereas the obtained p values were adjusted for multiple testing according to Benjamini and Hochberg (1995). The population from Port Edward (only N = 5) was excluded from this analysis. Owing to the clear differentiation between Halozetes populations, we calculated net average p distances between groups of sequences in MEGA v.7 (Kumar et al., 2016). Additionally, we calculated summary within-population statistics, such as genetic diversity indices (incl. number of segregating sites S, haplotype diversity H e , number of different haplotypes h, and nucleotide diversity π) in Arlequin v3.5.2.2.

Morphometric analyses
A total of 16 continuous variables were measured in 169 Halozetes capensis specimens from nine different spots (10 × 10 cm patch of algae or other substrate) originating from four sample locations (specimens of one sample location were pooled for the analyses). A total of 15 variables ( Fig. 2) were measured in 53 Fortuynia elamellata micromorpha from five different locations and in 62 Schusteria ugraseni specimens from four different localities.
Females and males were analyzed separately because of a pronounced, to a large part size-dependent sexual dimorphism, which is common in littoral oribatid mites and had also been revealed in preliminary analyses in the present study. Minimum, maximum, mean, standard deviation, and Data was size-corrected as described in Pfingstl et al. (2017), and multivariate analyses were performed on ln (x + 1) transformed raw as well as size-corrected data. If not indicated otherwise, analyses were performed with PAST 3.11 (Hammer et al., 2001).
Non-metric multidimensional scaling (NMDS) was used to visualize the pattern of variation within the analyzed species deriving from different sample locations. Linear discriminant analysis (LDA) was conducted in order to evaluate which variables were contributing most to the differences between the specimens from different locations and whether the same variables were important in both sexes. The performance of classification by LDA was tested by calculating the number of specimens correctly classified by leave-one-out cross-validation (LOO-CV) LDA, and the equality of means of the populations was tested by permutational multivariate analysis of variance (PERMANOVA). The differences in dispersion of specimens from each location were evaluated using the functions betadisper and permutest in the R package vegan (Oksanen et al., 2019).

Fortuynia e. micromorpha
TCS haplotype network analysis of COI sequence data of F. e. micromorpha populations showed relatively high genetic

Fortuynia e. micromorpha
The populations of F. e. micromorpha were overlapping in NMDS conducted on females as well as on males (Fig. 5). This was true for the raw data as well as for the sizecorrected data, although in size-corrected data of females a separation between the population from Umdloti and the populations from Umkomaas and Mtwalume (and partly also Sheffield) was present.  LDAs on both raw and size-corrected data of F. e. micromorpha females revealed a separation caused by the first two axes between the populations from Umdloti, Umkomaas, Mtwalume, and Sheffield, whereas the latter two formed a cluster in the size-corrected data (Fig. 5). The variables with highest loadings (and thus most responsible for separation between populations) in the raw data were nw da , nw dp , and dcg on axis 1 and ll, nw da , gl, al, and aw on axis 2. In the size-corrected data, highest loadings were present in nw da on axis 1 and in nw da and nw dp on axis 2 (Supporting Table S2). The power of classification by LDA was low: 10% of specimens in raw and 15% in size-corrected data were correctly classified. In accordance with this result, PER-MANOVA on both raw and size-corrected data revealed no significant differences between the female populations.
In contrast to the females, the males from the four populations of F. e. micromorpha were not clearly separated by LDA. In both raw and size-corrected data, small overlapping areas were present between the populations. In the raw data, variables with highest loadings were ll and dcg on axis 1 and nw dp and nw dm on axis 2. Variables with highest loadings in size-corrected data were ll and dcg on axis 1 and ll and nw dp on axis 2 (Supporting Table S2). The power of classification by LDA was better than in females; 39.39% were correctly classified in raw and 36.36% in size-corrected data. PER-MANOVA showed that there were significant differences (p < 0.01) between at least one of the male populations and the others in both raw and size-corrected data. Pairwise comparisons of the populations revealed significant differences (Bonferroni corrected p value p < 0.05) only between the populations from Umdloti and Umkomaas in size-corrected data.
Significant differences between the dispersion of populations could not be detected in any sex, or in raw or sizecorrected data.

Halozetes capensis
In NMDS on H. capensis, males overlapped in raw as well as in size-corrected data, but there was some separation between the females of different populations: the populations from Nature's Valley and Wilderness were largely separated from the populations from De Hoop and Kayser's Beach, and this separation was clearer in the size-corrected data (Fig. 6).
LDAs on females and males, based on raw as well as on size-corrected data, revealed similar patterns (Fig. 6): in all cases, the two populations from Nature's Valley and Wilderness were separated from the populations from De Hoop and Kayser's Beach on axis 1. Axis 2 separated the population from Kayser's beach from the other populations, whereas a small overlapping area remained in the males. The power of classification by LDA was also similar in all analyses: LDA on females correctly classified 67.8% in raw data and 61.02% in size-corrected data; in males, the percentages were higher with 72.73% correctly classified in raw data and 70.91% in size-corrected data. Thus, in both sexes, more specimens were correctly classified in the raw data.
Variables with highest loadings in females were db and nw dp for axis 1 and dga, al, and aw for axis 2 in the raw data, and gl, db, and nw da for axis 1 and dga and nw dp for axis 2 in the size-corrected data (Supporting Table S3). In males, the variables with highest loadings were db, nw dm , nw dp , and ddis for axis 1 and nw dp , dga, aw, and nw dm for axis 2 in the raw data. In size-corrected data, variables with highest loadings were db and nw da for axis 1 and nw dm and nw dp for axis 2 (Supporting Table S3).
PERMANOVA on both sexes, based on raw as well as on size-corrected data, showed that there were highly significant differences (p < 0.001) between at least one of the populations and the others. In pairwise comparisons between the female populations, significant differences (p > 0.05) were encountered between Nature's Valley and De Hoop and Nature's Valley and Kayser's Beach in both raw and size-corrected data. Pairwise comparisons of the male populations (raw as well as size-corrected data) showed significant differences (p < 0.01) between all possible pairings except between Wilderness and Nature's Valley, where no significant differences were present.
Significant differences in the dispersion of populations (p < 0.05) were found in the female populations between Wilderness and Nature's Valley and between Wilderness and De Hoop in the raw data. In size-corrected data of females, significant differences (p < 0.01) in dispersion were present only between Wilderness and Nature's Valley. In raw data of the male populations, the dispersion of the populations from Wilderness and Nature's Valley also differed significantly (p < 0.05); there were no significant differences in size-corrected data of males. In all mentioned cases, the population from Wilderness showed smaller dispersion than the other populations.

Schusteria ugraseni
The populations of S. ugraseni also overlapped to a large degree in NMDS on both sexes (Fig. 7). In the size-corrected data on females, the population from Winterstrand was separated from the other populations. The same trend was recognizable in raw data of females and in raw and size-corrected data of males, although the separation was less pronounced and the overlaps were larger in these cases.
LDAs on raw and size-corrected data on both sexes of S. ugraseni showed a clear separation of the population from Winterstrand and the other populations, which was furthermore more pronounced in males than in females (Fig. 7). The population from Winterstrand was in all LDAs separated from the others on axis 1. The power of classification by LDA was higher in females than in males, and in both sexes, it was higher in the size-corrected data. In females, 40% were correctly classified in raw data and 51.43% in sizecorrected data; in males, LDA correctly classified 14.81% in raw data and 25.93% in size-corrected data.
In females, the populations from Sheffield and Umkomaas clustered together in both raw and size-corrected data, albeit with only small overlapping areas (Fig. 7). They were separated from the population from Umdloti on axis 2. PER-MANOVA on raw and size-corrected data of female populations revealed significant differences (p < 0.01) between at least one of the populations and the others. In pairwise comparisons, there were significant differences (p < 0.05) between Winterstrand and Umdloti in raw data and between Winterstrand and each other population in size-corrected data.
In males, Sheffield and Umdloti clustered together in both raw and size-corrected data. The population from Umkomaas was separated from the other populations on axis 2, and the separation was more pronounced in the LDA on size-corrected data. Variables with highest loading in females were db i , nw dp , and gw for axis 1 and db i and nw dp for axis 2 in the raw data, and db i , cw, and gw for axis 1 and db i and nw dp for axis 2 in the size-corrected data. In males, the variables with highest loadings were bl, nw dp , and gw for axis 1 and gl, gw, al, and aw for axis 2 in the raw data. In size-corrected data, variables with highest loadings were dPtI and db i for axis 1 and nw dp and gl for axis 2 (Supporting Table S4). PERMANOVA on male populations revealed no significant differences in raw data and significant (p < 0.05) differences between at least one of the populations and the others in size-corrected data.
Significant differences (p < 0.05) in dispersion were present between the populations from Umkomaas and Winterstrand in both sexes, but only in the size-corrected data. In both sexes, the dispersion was larger in Winterstrand than in Umkomaas.

Genetic structures and their causes
Haplotype network analyses of COI sequence data of South African intertidal oribatid mite populations reveal clear but contrasting patterns. Populations of Fortuynia e. micromorpha show no geographic structuring and a single haplotype is shared by numerous individuals from all different locations, which indicates high levels of gene flow and connectivity across sampling localities. Populations of Halozetes capensis, on the other hand, exhibit a distinct geographic structure with no shared haplotypes, which points to restricted gene flow and isolation between the geographic locations. Paleoenvironmental processes may have played a considerable role in shaping this contrasting genetic pattern. During the last glacial maximum (LGM), ca. 22000 years ago, the Agulhas Bank became exposed due to a drop of 120 m in sea level and water temperatures decreased significantly (Romero et al., 2003). The Agulhas current slowed down and may have even ceased flow during winter (Dingle & Rogers, 1972;Hutson, 1980). Fortuynia e. micromorpha is a warm-adapted species (Pfingstl et al., 2021a) and the colder climate and weakening of the Agulhas current during the LGM surely shifted its distribution range to the north, whereas the cold-adapted H. capensis possibly expanded its occurrence to the east coast of South Africa. Consequently, South African F. e. micromorpha populations may have experienced a bottleneck during this period but expanded their distribution again with the following global warming and the strengthening of the Agulhas current. The star-like COI haplotype network supports the hypothesis of a recent population expansion for F. e. micromorpha. The cold-adapted H. capensis, on the other hand, was most likely pushed back to southern coastlines by the warmer climate following the LGM, but southern and southeastern populations experienced no dramatic temperature or ocean current changes and thus persisted throughout both periods. Both scenarios-i.e., demographic stable populations throughout the LGM and post-LGM recolonizationhave been suggested to shape populations structure in several other South African marine taxa (e.g., Heyden et al., 2007;Wood et al., 2017).
Although paleoenvironmental scenarios may explain regional persistence and recent expansion, they do not explain the significant differences in gene flow between Halozetes and Fortuynia. To answer that question, it is necessary to consider possible ways of dispersal and the factors influencing the dispersal potential. Rocky intertidal mites are air-breathing and mobile organisms that could possibly disperse via migration over land, but large dune fields and river mouths represent considerable barriers for large-scale population connectivity. From west to east along the South African shoreline, the intertidal habitat becomes progressively less rocky and is interspersed with numerous sandy beaches (Heyden et al., 2008). Hence, theoretically, the Western Cape populations of H. capensis experience fewer overland barriers than their eastern counterparts, F. e. micromorpha. Nonetheless, the exact opposite is reflected by the present genetic data, which largely contradicts successful over land migration. Moreover, given the small size of the animals, it seems unlikely that they actively migrate for several kilometers from one population to the other with the above-mentioned barriers in between. Our genetic data also indicate that migration even on a small local scale is rather uncommon as certain samples from single spots on a location only showed a single haplotype (e.g., ZA_64, JA_05; see supporting Table S5). Nevertheless, in other locations, as for example De Hoop, we found two different haplotypes being shared among four different sample spots at this location (see Supporting Table S5), and these spots were placed within in a radius of approx. 100 m. This shows that genetic exchange between the spots of a single location can take place and small-scale migration could still be responsible.
Other possible ways for overland dispersal for intertidal oribatid mites include passive wind drifting or transport in the feathers of shorebirds. Even though wind dispersal of a few terrestrial oribatid mites has been observed, the majority of these mites consisted of species usually living in the canopy of trees (Lehmitz et al., 2011). Intertidal oribatid mites, on the other hand, dwell in algae growing in crevices of rocks and are thus more or less sheltered from winds and less prone to be blown away. Furthermore, they have not been found in the plumage of shorebirds yet and these mites do not show any morphological adaptations to attach to a bird's body. Consequently, neither wind nor bird dispersal could explain the contrasting population patterns observed in this study.
Hence, a more likely dispersal mode for the investigated species might be passive drift on ocean currents. Pfingstl et al. (2021a) demonstrated that distribution ranges of South African intertidal oribatid mites coincide with marine biogeographic regions and that oceanic climate influences occurrence patterns. Several authors (e.g., Pfingstl, 2013b;Schatz, 1991) argued that drifting along ocean currents is the main mode of long-distance transport for intertidal mites and it was shown that some Caribbean species would theoretically be able to survive more than a month drifting in the Gulf Stream (Pfingstl, 2013b). This could also be true for the southern African coasts, where the population connectivity of marine species is strongly linked to oceanic currents. In this context, four major gene flow scenarios have been identified: (I) strong northward exchange with the Benguela Current on the west coast, (II) strong southward flow with the Agulhas Current on the east coast, (III) some bidirectional connectivity inshore of the Agulhas Current on the southeast coast, and (IV) bidirectional gene flow on the south coast (Teske et al., 2011). The lack of genetic differentiation in F. e. micromorpha populations could thus be related to the Agulhas current on the east coast. The Agulhas current flows close inshore with high velocity off northern KwaZulu-Natal (Schumann, 1987) and could easily transport mites that were washed off the rocks to another location on the coast. Stochastic but frequent transport along this strong ocean current could explain the low level of genetic structuring found in the populations of F. e. micromorpha. Unfortunately, the inference about the direction of gene flow (e.g., by means of an isolation model with migration) was not possible, due to the low resolution of the selected markers in this study. The investigated H. capensis populations, on the other hand, occur in the warmtemperate Agulhas Ecoregion, where the Agulhas current moves well offshore following the edge of the Agulhas Bank (Grundlingh, 1983). The deflection of the Agulhas current is supposed to be responsible for the genetic discontinuity found in certain shrimp (Wood et al., 2017) and it could also be responsible for the lack of gene flow between Halozetes populations. Nonetheless, a South African bluntnose klipfish species is supposed to utilize smaller counter-currents that develop inshore of the Agulhas current for dispersal along the southeast coast (Heyden et al., 2008). In contrast to fish, intertidal oribatid mites do not possess highly dispersive planktonic larvae. Additionally, to properly migrate over sea, they need to be washed out far enough to reach large oceanic circulation systems. Since this happens rather accidentally than actively (e.g., during heavy storms), local isolation is maintained to a certain extent.
We observed that populations of H. capensis from Wilderness and Nature's Valley are closer related to the population from Kayser's Beach than to the populations from De Hoop. De Hoop lies approx. 200 km to the west of Wilderness and Nature's Valley, whereas Kayser's Beach is located approx. 400 km to the east of these two locations. The Agulhas current is still close to the shore in the area of Kayser's beach and therefore westward transport of mites to the area near Port Elizabeth, where it finally moves offshore, could still be possible. This would explain why phylogenetic relatedness in H. capensis does not coincide with geographic distance.
Two factors may additionally contribute to the diverging genetic patterns observed in this study. Firstly, the Caribbean Fortuynia atlantica Krisper & Schuster, 2008 was observed to show a special floating behavior when suddenly washed off the substrate (Pfingstl, 2013a), perhaps a similar behavior is present in F. e. micromorpha, which would clearly facilitate hydrochorous dispersal and enhance population connectivity. Secondly, the genus Halozetes comprises intertidal and typical terrestrial species (Marshall & Convey, 2004) indicating an evolutionary weaker bond to the intertidal environment. Consequently, Halozetes may be less resistant to long-term submergence and thus less prone to dispersal via ocean currents. Nevertheless, experimental studies investigating the behavior of F. e. micromorpha and the tolerance of H. capensis to submergence in salt water are necessary to confirm these hypotheses.

Morphological variation
Morphometric investigations revealed considerable morphological homogeneity in the studied F. e. micromorpha populations. Variation does not reflect a geographic pattern and is consistent with molecular genetic results indicating that the found homogeneity is most likely a result of the high level of gene flow and connectivity between the populations. Similar results are shown for Schusteria ugraseni populations, but the population from Winterstrand is clearly separated from the others in the morphometric analyses. Winterstrand is located more than 320 km southwest from the other sample locations whereas the distances between the latter only range from 25 to 60 km. Therefore, the morphological differentiation of the population from Winterstrand is most likely a result of geographic separation and accompanying restricted gene flow. Unfortunately, no genetic data for this species is available allowing further interpretation and comparison with the other taxa.
Despite large overlaps, two morphological groups are present among H. capensis populations, one group comprises Nature's Valley and Wilderness and the other consists of the De Hoop and the Kayser's Beach populations. The first group clearly reflects geography and genetic relatedness; the latter group, however, does not. De Hoop and Kayser's Beach are more than 600 km apart and the population from Kayser's Beach is most closely related to the populations from Wilderness and Nature's Valley. The morphological similarity between De Hoop and Kayser's Beach can therefore not be explained by geographic proximity or genetic relatedness. Pfingstl et al. (2021b) suggested that morphological variation contrasting with genetic data in intertidal oribatid mites could be a result of phenotypic plasticity caused by diverging local environmental factors. Accordingly, it is possible that similar environmental properties of the De Hoop and the Kayser's Beach sample sites resulted in similar phenotypes of the populations. Further on-site studies are necessary to identify these possible similar environmental factors and to verify the suggested correlation.