Shedding light on the Ophel biome: the trans-Tethyan phylogeography of the sulfide shrimp Tethysbaena (Peracarida: Thermosbaenacea) in the Levant

Introduction

Groundwater fauna (stygofauna) is characterized by short-range endemism and high species crypticity (Eme et al., 2013; Zagmajster et al., 2014). The unique suite of troglomorphic traits (e.g., loss of pigment, reduced eyes) characterizing stygobionts often hinders distributional studies due to the highly convergent characteristics that can obscure taxonomic relationships (Juan et al., 2010; Porter, 2007). As a result, molecular phylogenetic tools have been extensively used over the last two decades to infer stygofauna biogeographies and the underlying processes shaping them (e.g., Abrams et al., 2019; Asmyhr et al., 2014; Bauzà-Ribot et al., 2012; Bradford et al., 2010; Cánovas et al., 2016; Cooper et al., 2023; Finston et al., 2004; Guy-Haim et al., 2018; Jaume, 2008; Jurado-Rivera et al., 2017; Marin, Krylenko & Palatov, 2021; Matthews et al., 2020).

Thermosbaenacea is a small order of peracarid crustaceans comprising 36 unique and highly specialized species adapted to extreme aquatic environments, including spring-fed subterranean lakes and thermal springs, with their core populations found deep underground in the inaccessible phreatic waters (Por, 2014; Wagner, 1994; WoRMS, 2023). Anoxic, sulfide-rich environments are favorable to Thermosbaenacea—often feeding on bacterial mats formed by sulfide-oxidizing bacteria—thus termed “sulfide shrimp” by Por (2014). Based on their distribution, it was assumed that the ancestral habitat of the thermosbaenaceans is the Tethys—an ancient ocean that was situated between Laurasia to the north and Gondwana to the south and existed from the late Paleozoic to the early Cenozoic eras (ca. 250–5 MYA). Nowadays, Thermosbaenacea representatives, largely relict and geographically isolated species, are found around the Mediterranean, the Arabian Sea, the Indian Ocean, Mid-East Atlantic, and the Caribbean Sea (Hou & Li, 2018; Wagner, 1994). Among thermosbaenaceans, Tethysbaena (family: Monodellidae) is the most speciose and widespread genus, comprising 27 species in seven species-groups (Wagner, 1994; Wagner & Bou, 2021). Only a few of the Tethysbaena species-groups were analyzed and supported by molecular phylogenetic tools (Cánovas et al., 2016; Wagner & Chevaldonné, 2020).

Two species of Tethysbaena are known from Israel: T. relicta (Por, 1962) (formerly Monodella relicta) and T. ophelicola (Wagner, 2012). Initially, fragments of T. relicta were found in the hot sulfide-rich spring Hamei Zohar by the Dead Sea in Israel (Por, 1962). Later, a few specimens of the same species were occasionally found in the thermohaline spring En-Nur, on Lake Kinneret shore, a few hundred kilometers to the north (Dimentman & Por, 1991), thus inferring that T. relicta inhabits the whole groundwater system of the Dead Sea-Jordan Rift Valley aquifer. T. ophelicola was found in the karstic underground basin near Ramla, named Ayyalon-Nesher-Ramla complex (Por, 2014; Por et al., 2013; Wagner, 2012), 60 km west of the Dead Sea-Jordan Rift Valley, beyond the water divide of Israel, thus belonging to a different watershed system. Similarly to other Thermosbaenacea species, the Levantine Tethysbaena inhabit environments characterized by high temperatures and sulfide-rich waters—31 °C in Hamei-Zohar (Por, 1962), 30–31 °C with 60–190 ppm sulfate in En-Nur (Beeri-Shlevin et al., 2022; Tsurnamal, 1978), and 28.8–29.6 °C with 160–305 ppm sulfate concentrations in the Ayyalon cave complex (Langford et al., 2022; Por et al., 2013).

Based on synapomorphies of the antennular inner flagellum and maxilliped macrosetae (Wagner, 1994), it was hypothesized that together with other closely allied species (one species from Somalia (Chelazzi & Messana, 1982), four species from Oman (Wagner, 2020), one species from Yemen (Wagner & Van Damme 2021)), T. relicta and T. ophelicola form the “T. relicta-group” (Wagner, 2012), suggesting a recent common ancestor. An alternate hypothesis can be drawn from the phylogenetic analysis of the prawn Typhlocaris (Guy-Haim et al., 2018), preying on Tethysbaena in Ayyalon and En-Nur (Tsurnamal, 1978; Tsurnamal, 2008; Tsurnamal & Por, 1971; Wagner, 2012). Four Typhlocaris species are known, two of which co-occur with Tethysbaena in Israel: Ty. galilea inhabiting En-Nur spring (Calman, 1909; Tsurnamal, 1978) and Ty. ayyaloni from the Ayyalon cave (Tsurnamal, 2008). The two additional Typhlocaris species are Ty. salentina from Southeastern Italy (Caroli, 1923; Froglia & Ungaro, 2001) and Ty. lethaea from Libya (Parisi, 1921). The molecular phylogeny of Typhlocaris species showed that Ty. ayyaloni (Israel) and Ty. salentina (Italy) are more closely related to each other than either of them is to Ty. galilea (Israel) (Guy-Haim et al., 2018). Accordingly, we can hypothesize a similar phylogeographic pattern of the Levantine Tethysbaena, where T. ophelicola would be more closely related to the Mediterranean species (“T. argentarii-group”) than to T. relicta.

Along the biospeleological research of the Thermosbaenacea and other phyla of subterranean crustacean fauna represented in the Dead Sea Rift Valley (Syncarida, the families Bogidiellidae, Typhlocarididae, Niphargidae, Sphaeromatidae, and Cirolanidae), three paradigms have been proposed to explain their origins. Accordingly, Tethysbaena colonization could be attributed to: (1) Pliocene pre-glacial (Piacenzian, 3.600–2.588 MYA) marine transgression (Por, 1963), penetrating the Dead Sea Rift Valley from the Mediterranean; (2) mid-Miocene Tethys transgression into the Levant (16–10 MYA) (Dimentman & Por, 1991; Por, 1987); and (3) The Ophel Paradigm that offered a conceptual framework, within which these styobionts are inhabitants of the ancient chemosynthetic Ophel biome, dating back at least to the Cambrian (>485 MYA) (Por, 2011).

Using a molecular clock approach, Guy-Haim et al. (2018) estimated the divergence time of the Typhlocaris species. They based their analysis on a calibration node inferred from a regional geological event—the end of the marine connection between the Mediterranean Sea and the Dead Sea-Jordan Rift Valley, marked by the top of Bira formation, dated to 7 MYA (Rozenbaum et al., 2016), separating Ty. galilea and the Typhlocaris ancestor. The inferred divergence time of Ty. ayyaloni and Ty. salentina was 5.7 (4.4–6.9) MYA, at the time of the Messinian Salinity Crisis (5.96–5.33 MYA), when the Mediterranean Sea desiccated (Garcia-Castellanos & Villaseñor, 2011; Gargani & Rigollet, 2007) and lost almost all its Miocene tropical fauna (Por, 1987; Por, 1989). It is therefore an open question as to whether the same vicariant events have shaped the biogeographies of both predator (Typhlocaris) and prey (Tethysbaena) subterranean crustaceans.

The main objectives of our study were to (1) reveal the phylogenetic relatedness of the Levantine Tethysbaena species, and use these patterns to (2) infer the geological and evolutionary processes that have shaped their divergence patterns. To do so, we used the mitochondrial COI and nuclear 28S rRNA genes to construct the phylogeny of T. relicta and T. ophelicola and estimated their divergence ages by applying molecular clock approach.

Materials & Methods

Sampling sites, specimen collection and identification

Specimens of T. ophelicola were collected by a hand pump from the surface water of the inner pool (maximal depth ca. 2 m) of the Levana cave (31.9223°N, 34.8942°E), part of the Ayyalon-Nesher-Ramla complex (Fet et al., 2017) (Fig. 1).

Specimens of T. relicta were collected by a hand pump from the overlying water near the bottom of an artificial tunnel (maximal depth ca. 0.6 m) near the Dead Sea Shore penetrating the Judea Group aquifer, 6.5 km north of Hamei-Zohar (31.2232°N, 35.3547°E) (Fig. 1). The locus typicus of T. relicta, the thermal spring of Hamei-Zohar (Por, 1962), is no longer accessible since the 1970s, as hotels were built on the spring area. Nonetheless, since Dimentman & Por (1991) found scattered T. relicta specimens in En-Nur, 150 km north of the type locality Hamei-Zohar, we assumed that the new Tethysbaena record, located along the Dead Sea Jordan Rift, is of T. relicta, and validated it using morphological characterization.

Part of the collected specimens was preserved in 70% ethanol and the other in absolute ethanol for morphological and molecular analyses, respectively. Species identification of T. ophelicola and T. relicta was performed using a stereomicroscope (SZX16, Olympus, Japan) following the identification keys in Por (1962), Wagner (1994) and Wagner (2012).

Results

Morphological identification

Specimens of T. relicta collected from the Dead Sea tunnel were similar to the specimens from the type locality, Hamei-Zohar thermal spring, described by Por (1962), and included males, with no ovigerous or brooding females (Fig. 2A). The average length (excluding antennae) was 2,104 ± 181 µm (n = 5, ±SD, hereafter). The following morphological features characterized the specimens as belonging to T. relicta: eight segments in the main flagellum (endopodite) of antenna 1; seven terminal plumidenticulate macrosetae were present on the maxilliped; the uropod included five medial plumose macrosetae, 11–13 plumose macrosetae in the endopodite, and 16–19 macrosetae in the second segment of the exopodite. The mean width:length ratio of the telson was 1.15.

Specimens of T. ophelicola from Levana cave were similar to the specimens from Ayyalon cave described by Wagner (2012), and included males, brooding females and postmarsupial juveniles (Figs. 2B–2D). The average length (excluding antennae) was 2,276 ± 380 µm in males (n = 5) and 2,620 ± 139 µm in females (n = 5). The following morphological features were found: seven segments in the main flagellum (endopodite) of antenna 1; seven terminal plumidenticulate macrosetae were present on the maxilliped; uropod included 4 medial plumose macrosetae and 18–22 plumose macrosetae in the endopodite and the second segment of the exopodite. The mean width:length ratio of the telson was 1.10.

Molecular phylogenetic analysis

The DNA barcode consisting of a fragment of 624–641 bp of the COI gene was sequenced from 6 specimens of T. ayyaloni and T. relicta. Sequences were deposited in NCBI GenBank under accession numbers OR189199 –OR189204. The phylogenetic analysis included 16 additional Tethysbaena sequences and one Halosbaena tulki sequence as an outgroup (Table 1). The overall COI alignment was 691 bp long, with 227 parsimonious informative sites.

The amplified 28S rRNA gene sequences of T. ayyaloni and T. relicta, deposited in NCBI GenBank under accession numbers OR790561 –OR790565, were 834–885 bp long. The phylogenetic analysis included an additional Tethysbaena sequence (T. argentarii) and three Halosbaena sequences as an outgroup (Table 2). The overall 28S rRNA alignment was 908 bp long, with 498 parsimonious informative sites.

Table 2:

28S rRNA sequences of Tethysbaena and outgroup included in the phylogenetic analysis.

	Species	Accession number	Locality	Reference
1	Tethysbaena relicta	OR790563.1	Dead Sea tunnel, Israel	This study
2	Tethysbaena relicta	OR790564.1	Dead Sea tunnel, Israel	This study
3	Tethysbaena relicta	OR790565.1	Dead Sea tunnel, Israel	This study
4	Tethysbaena ophelicola	OR790561.1	Levana cave, Israel	This study
5	Tethysbaena ophelicola	OR790562.1	Levana cave, Israel	This study
6	Tethysbaena argentarii	DQ470654.1	Italy	Stenderup, Olesen & Glenner (2006)
7	Halosbaena fortunata (outgroup)	KT984015.1	Spain	Page et al. (2016)
8	Halosbaena daitoensis (outgroup)	KT984014.1	Japan	Page et al. (2016)
9	Halosbaena tulki (outgroup)	KT984043.1	Australia	Page et al. (2016)

DOI: 10.7717/peerj.16690/table-2

The COI-based ML and Bayesian phylogenetic analyses showed similar tree topologies (Fig. 3). The Levantine Tethysbaena species from Israel present polyphyly, where T. ayyaloni is positioned within a Mediterranean clade (including T. scabra from the Balearic Islands, T. ledoyeri from Southern France and T. argentarii from Italy) with 100% bootstrapping support and 0.99 posterior probability, and T. relicta clusters with Tethysbaena sp. from Oman (100% bootstrapping support and 1.00 posterior probability), and the Dominican Republic (87%/0.83 bootstrapping support/posterior probability), forming the Arabian-Caribbean clade. The Atlantic T. atlantomaroccana is a sister taxon to the Mediterranean clade species, although with a lower support/probability. The other Moroccan Tethysbaena species from Tasla and Lamkedmya were in a more basal position but showed lower bootstrapping support (<50%).

The 28S rRNA-based ML and Bayesian phylogenetic analyses (Fig. 4) included only a subset of the species in the COI phylogenetic analysis. Similar to the COI-based phylogenetic analysis, it showed that the Levantine T. ayyaloni is a sister taxon to the Mediterranean clade species T. argentarii from Italy (86%/1.00 bootstrapping support/posterior probability), whereas T. relicta is more distant (100%/1.00 bootstrapping support/posterior probability).

Divergence time estimation

Effective sample size (ESS) values were at least 436 and 356 for the posterior and prior statistics, respectively, 1,738 for the likelihood statistic, and greater than 1,400 for all MRCA times estimates, suggesting good mixing and an effective MCMC sampling of the posterior distribution.

We estimated the ages for eight nodes (Table 3, Fig. 5). The youngest node was the most recent common ancestor of T. leyoderi from Southern France and T. scabra from the Balearic Islands, which returned a mean estimate at 8.31 MYA. The next mean estimate is the divergence of T. ophelicola from the clade of T. leyoderi and T. scabra, dated to 9.46 MYA. The mean age of the most common ancestor of all Mediterranean Tethysbaena was 10.71 MYA. The most recent ancestor of the Mediterranean clade and T. atlantomaroccana from Morocco was dated to 32.41 MYA. The mean age of the node linking T. relicta from the Dead Sea-Jordan Rift Valley and Tethysbaena from Oman was 20.13 MYA. The node of the most recent common ancestor of T. relicta, Tethysbaena from Oman, and the Tethysbaena from the Dominican Republic had a mean estimate of 35.84 MYA. The mean age for the node linking the Arabian-Caribbean clade (T. relicta + Tethysbaena sp. Oman + Tethysbaena sp. Dominican Republic) with the Mediterranean-Atlantic clade (T. scabra + T. leyoderi + T. ophelicola + T. argentarii + T. atlantomaroccana) was 40.42 MYA. The estimate for the root node linking Tethysbaena and Halosbaena was 79.96 MYA.

Table 3:

Divergence times for Tethybaena species as estimated by the Bayesian evolutionary analysis method calculated using the COI gene molecular evolution based on Cánovas et al. (2016) and Bauzà-Ribot et al. (2012).

Node ages and highest posterior density ( ±95% HPD) ranges are given in million years round.

	Clade divergence (nodes)	Node age (MYA) (95% HPD range)	Geological period
1	T. scabra—T. leyoderi	8.31 (10.15–3.97)	Miocene
2	T. scabra + T. leyoderi—T. ophelicola	9.46 (14.20–5.71)	Miocene
3	T. scabra + T. leyoderi + T. ophelicola—T. argentarii	10.71 (16.27–6.04)	Miocene
4	T. scabra + T. leyoderi + T. ophelicola + T. argentarii— T. atlantomaroccana	32.41 (47.53–18.37)	Oligocene
5	T. relicta—Tethysbaena sp. (Oman)	20.13 (41.69–13.25)	Miocene
6	T. relicta + Tethysbaena sp. (Oman)—Tethysbaena sp. (Dominican Republic)	35.83 (51.41–22.16)	Eocene
7	T. scabra + T. leyoderi + T. ophelicola + T. argentarii + T. atlantomaroccana—T. relicta + Tethysbaena sp. (Oman) + Tethysbaena sp. (Dominican Republic)	40.42 (56.09–25.72)	Eocene
8	Tethysbaena—Halosbaena	79.96 (137.8–32.68)	Cretaceous

DOI: 10.7717/peerj.16690/table-3

Discussion

In his monography on Thermosbaenacea, Wagner (1994) divided the Monodellidae family to two genera, the monotypic Monodella and the speciose Tethysbaena, which he named after the ancient Tethys Sea and the Greek word “βαινɛιν” (meaning “to walk”), referring to these animals as “walkers of the Tethys Sea”. He noted that although there is a great similarity among the different species, six species-groups can be identified based on morphological characters. With the later finding of T. exigua from Southern France, a seventh group was established (Wagner & Bou, 2021). Here, we analyzed the phylogenetic relatedness and divergence times of the two Levantine Tethysbaena species found in Israel: T. relicta from the Dead Sea-Jordan Rift Valley, and T. ophelicola, from the Ayyalon-Nesher-Ramla cave complex in central Israel.

According to Wagner (2012) and Wagner & Van Damme (2021), both Levantine species belong to “T. relicta-group” (together with four species from Oman, one species from Somalia and one species from Yemen), implying that these are sister taxa sharing a most recent common ancestor. Our phylogenetic results, based on the mitochondrial COI and the nuclear 28S rRNA genes, reject the morphology-based cladistics. The COI-based phylogeny supports the hypothesis suggesting that T. relicta shared an ancestor with Tethysbaena species from Oman and Dominican Republic, whereas the circum-Mediterranean species (including T. ophelicola) share another ancestor. Indeed, discrepancies between morphological cladistics and molecular phylogeny are common in cave fauna and were often attributed to their convergent troglomorphic traits (Bishop & Iliffe, 2012; Juan et al., 2010; Liu et al., 2017; Porter, 2007; Sendra et al., 2021).

Three paradigms determining the origin of the Thermosbaenacea and other phyla of subterranean crustaceans represented in the Dead Sea-Jordan Rift Valley and around the Mediterranean were defined. The earlier paradigm suggested that the Levantine Tethysbaena, among other subterranean salt-water fauna, have resulted from a late Pliocenic pre-glacial (Piacenzian, 3.600–2.588 MYA) marine transgression (Fryer, 1964; Hubault, 1937; Por, 1963). A narrow gulf penetrated into the coastal line near the present-day mount Carmel and then bent southwards along the Dead Sea-Jordan Rift Valley reaching a basin that extended south of the present Dead Sea (Picard, 1943). According to this paradigm, the Pliocenic Mediterranean was still inhabited by a very large number of Tethys remnants, including thermosbaenaceans, that were stranded in the Rift Valley and around the Mediterranean.

Por (1986) rejected the first paradigm, noting that the Pliocenic Mediterranean no longer contained the tropical fauna, including the Tethysbaena ancestor, and that the short-lived Pliocene transgression did not establish viable marine environments. Instead, he posited that these species represent marine fauna colonized during a Miocenic transgression (16–10 MYA) (Bar & Zilberman, 2016; Buchbinder & Zilberman, 1997), the last time that tropical sea penetrated inland in the Levant, and left stranded following a late Miocene regression (6–5.3 MYA) (Dimentman & Por, 1991; Por, 1987; Por, 1989). This second paradigm was supported by Guy-Haim et al. (2018) who applied a molecular clock approach to estimate the divergence time of the co-occurring caridean prawn Typhlocaris. Their assessment was based on a calibration node inferred from the end of the marine connection between the Mediterranean Sea and the Dead Sea-Jordan Rift Valley, marked by the top of Bira formation, dated to 7 MYA (Rozenbaum et al., 2016). They inferred a divergence time of Typholocaris from Ayyalon cave and Italy of 5.7 (4.4–6.9) MYA, at the time of the Messinian Salinity Crisis. During this event, the African plate moved towards the Euro-Asian plate, closing the Straits of Gibraltar and temporarily isolating the Mediterranean Sea from the Atlantic Ocean (Garcia-Castellanos & Villaseñor, 2011; Krijgsman et al., 1999). As a result, the Mediterranean Sea partly desiccated and transformed into small hypersaline basins, losing almost all its Miocenic tropical fauna, including those able to colonize subterranean waters (Carnevale et al., 2019; Delić et al., 2020; Por, 1987; Por, 1989).

With the discovery of the Ayyalon cave system and its endemic stygofauna in 2006, a third paradigm known as “the Ophel Paradigm” was developed by Por (2007). He identified the “Ophel” as a continental subterranean biome, subsisting on chemoautotrophic bacterial food, independently of the exclusive allochthonous epigean food of photoautotrophic origin. Within this biome, Tethysbaena are primary consumers, presenting a typical feeding behavior of upside-down swimming-gathering of sulfur bacteria or bacterial mats (Por, 2011; Wagner, 2012). Following the development of the new chemosynthetic-based biome paradigm, Por presented an alternative to the Tethys stranding paradigm, stating that the “Ophel paradigm falsified first of all my own, previously held views” on the diversification of the subterranean fauna in the Levant (Por, 2011). He noted that the pre-Messinian fauna of the fossiliferous taxa of the foraminiferans, the mollusks and the teleost fishes was similar to the recent Red Sea fauna or different only at the species level, and there is no indication for extinction of crustaceans during the Tertiary, thus inferring that the origin of the subterranean Levantine fauna is of earlier origin (Por, 2010). Por suggested that the Ophelic biome is possibly at least as old as the Cambrian, which had a diverse aquatic crustacean and arthropodan palaeofauna, including Thermosbaenacea (Por, 2011).

Cánovas et al. (2016) assessed the divergence time of the Western Mediterranean Tethysbaena, T. scabra from the Balearic Islands and T. argentarii from Italy using the COI gene. They based the nucleotide substitution rates on the mean rate estimated for a co-occurring anchialine stygobiont amphipod Metacrangonyx longipes, 1.32% per lineage and million years (95% CI [0.89–1.95]) (Bauzà-Ribot et al., 2012) and estimated the divergence time of T. scabra and T. argentarii to the early Tortonian, 10.7 MYA. Following Cánovas et al. (2016), we have used the COI gene to assess the divergence times of the Levantine Tethysbaena, T. relicta and T. opehlicola, and additional Tethysbaena species from around the Mediterranean, Arabian, and Caribbean Sea, using the Australian Halosbaena as an outgroup. Although the evolutionary rate used for calculating the Tethysbaena divergence closely matches the typical arthropod protein-coding mitochondrial substitution rates (1.15%/MY, Brower, 1994), it should be noted that the results based on this rate estimation should be treated with caution, as Metacrangonyx and Tethysbaena are only distantly related, and rates of molecular evolution tend to vary across genes, lineages and timescales (Hipsley & Müller, 2014).

Our analysis shows that the divergence times of Tethysbaena species are earlier than those of Typhlocaris species, pre-dating the upper-Miocene Messinian Salinity Crisis. Most divergence events occurred in the Miocene and Oligocene. The Dead Sea-Jordan Rift Valley T. relicta shares a most recent common ancestor with Tethysbaena from the Arabian Sea (Oman), dated to the early Miocene, 20.13 MYA (with 95% HPD of 41.69–13.25), corresponding with the Oligo-Miocene rift-flank uplift of the Arabian plate during the formation of the Red Sea and Gulf of Aden (34-20 MYA) (Omar & Steckler, 1995; Stern & Johnson, 2010). Both T. relicta and the Tethysbaena from Oman separated from the Caribbean Tethysbaena during the Eocene-Oligocene transition (38-30 MYA), when global cooling and tectonic uplift caused sea level decline and led to the establishment of the modern Caribbean Seaway (Iturralde-Vinent & MacPhee, 1999; Iturralde-Vinent, 2006; Weaver et al., 2016).

The most recent common ancestor of the Mediterranean Tethysbaena species—T. ophelicola from the coastal plain of Israel, T. scabra from the Balearic Islands, T. ledoyeri from Southern France, and T. argentarii from Italy—dated to the Tortonian in the Mid Miocene, 10.71 MYA (with 95% HPD of 6.27–6.04) as was previously found by Cánovas et al. (2016). The Ayyalon cave Tethysbaena, T. ophelicola, separated from other Mesogean (emerging Mediterranean) species around that time, 9.46 MYA (with 95% HPD of 14.20–5.71). The thermal water of the Ayyalon cave complex is part of the Yarkon-Taninim aquifer (Weinberger et al., 1994). During Oligocene-Miocene regressions (28-6 MYA), canyons were entrenched along the Mediterranean Sea shoreline, serving as major outlets of the Yarkon-Taninim aquifer, potentially forming anchialine karst caves (Frumkin, Dimentman & Naaman, 2022; Laskow et al., 2011). The ancestral habitat of Thermosbaenacea was first suggested to be of shallow marine origin (Fryer, 1964), stranded in the interstitial or inland groundwaters during periods of marine regression (Por, 1963; Stock, 1976; Tsurnamal & Por, 1971). An alternative hypothesis was proposed by Iliffe et al. (1984), who suggested that these stygobionts originate in the deep sea. This was later criticized by Stock (1986), among other reasons, because Thermosbaenacea colonizing marine anchialine caves have no relatives in the deep sea, but exhibit closer affinities to littoral interstitial species. Page et al. (2016) hypothesized that the ancestral habitats of Thermosbaenacea are Tethyan anchialine caves. Following our findings, we assume that the ancestor of T. ophelicola inhabited coastal anchialine caves in the Miocenic Tethys.

The most recent common ancestor of the Mediterranean and the Arabian-Caribbean clades of Tethysbaena is dated to the upper Eocene (40.42 MYA). During that period, the collision between the Arabian Plate and the Eurasian Plate resulted in the uplift of the Zagros Mountains in Iran (Mouthereau, Lacombe & Vergés, 2012). These mountain ranges acted as barriers, further isolating the Arabian Sea from the Mediterranean region (Sanmartín, 2003). The oldest, root node (Tethysbaena-Halosbaena) dated to 79.96 MYA (with the caveat of a low posterior probability and a large 95% HPD interval, 137.8–32.68 MYA). Page et al. (2016) established the phylogeny and divergence dates of the thermosbaeancean Halosbaena. They used the Tethysbaena-Halosbaena divergence as a calibration node, based on the presence of a continuous band of ocean crust through the length of the North Atlantic, indicating the maximum extent of the Tethys and the final opening of the Atlantic, dated to 107.5 MYA (with 95% HPD of 125–90). Thus, Tethysbaena ancestor in both our analysis and in Page et al. (2016) dates to the Cretaceous. The validity of the Paleozoic Ophel-driven hypothesis is also undermined by the deep phylogeny of peracaridean orders based on the small-subunit (SSU) rRNA gene, which showed that the thermosbaenacean lineage does not occupy a basal position relative to other peracarids (Spears et al., 2005). Numerous terrestrial sites having rich sulfidic waters have been discovered in the last century, and also recently (Nicolosi et al., 2022; Popa et al., 2019), after Por (2007); Por (2011) presented his “Ophel paradigm”. In general, we don’t see much evidence for a terrestrial Ophelic biome as old as the Cambrian. Rather, biogeographic, molecular, and hydrogeologic studies indicate that Cenozoic stranding events seem to be responsible for much of the observed subterranean fauna in these habitats (e.g., Delić et al., 2020; Esmaeili-Rineh et al., 2015; Frumkin, Dimentman & Naaman, 2022; Guy-Haim et al., 2018; Page et al., 2016).

Overall, the molecular clock-based divergence patterns presented here do not support the previously proposed hypotheses regarding the origins of the Levantine Tethysbaena. Instead, we infer a complex, two-stage colonization pattern of the Tethysbaena species in the Levant: (1) a late Oligocene transgression event, through a marine gulf extending from the Arabian Sea in the East to the Sea of Galilea in the west, leading to the colonization of T. relicta in the Dead Sea-Jordan Rift Valley, and (2) T. ophelicola, originating from the Mesogean Sea ancestor, inhabited anchialine caves in the coastal plain of Israel during the Mid-Miocene. Our results also show that the Cretaceous Tethysbaena ancestor first established in present-day Morocco, and then diverged into two groups. The first is a Tethyan group including Oman, the Dead Sea-Jordan Rift Valley and the Caribbean Sea. The second group formed around the emerging Mediterranean Sea, in its marginal aquifers, including Ayyalon, Southern France, Italy and the Balearic Islands.

Conclusions

Our multi-marker phylogenetic results reject the morphology-based cladistics and suggest that T. relicta shared a most recent common ancestor with Tethysbaena species from Oman and Dominican Republic, whereas the circum-Mediterranean species, including T. ophelicola, shared another ancestor. The molecular dating analysis suggests a two-stage colonization of the Tethysbaena species in the Levant, explaining their distant origins: a late Oligocene transgression leading to the colonization of T. relicta in the Dead Sea-Jordan Rift Valley, and a Miocene transgression in the Mediterranean region followed by a marine regression, stranding T. ophelicola in the coastal plain of Israel. The speciose Tethysbaena provides an exquisite opportunity for testing paleogeographic paradigms. Here we analyzed the phylogenetic relationships and divergence of nine out of twenty-seven known Tethysbaena species using the mitochondrial COI gene and the nuclear 28S rRNA gene. Future studies should examine additional species utilizing more genes or complete genomes to further unveil the phylogeny and biogeography of this unique and understudied group of ancient subterranean crustaceans.

The study of these subterranean species is not only an opportunity to broaden our understanding of paleogeography. It is also paramount for the protection of the hidden biodiversity found in these largely inaccessible habitats, which is nonetheless increasingly affected by human activity. Extraction of groundwater for irrigation and other uses, pollution, as well as quarrying, mining, and above-ground development may put these underground ecosystems at severe risk. The unique and often endemic nature of stygobiont species makes them even more prone to extinction, and extensive exploration of this under-explored biome, worldwide, is necessary in order to gain understanding and appreciation of the hidden biodiversity underground –an understanding that may pave the way for conservation of these species and their ecosystems.

Supplemental Information