Introduction

The present-day species- and genus-level diversity on earth is mainly the result of the evolution during Cenozoic, when biodiversity recovered from the Cretaceous-Tertiary mass extinction event, which eradicated three-quarters of all animal and plant species. During the Cenozoic new lifeforms flourished; often from lineages that had played only a minor role before the mass extinction; like mammals, birds and teleost fishes (Feduccia, 1995; Friedman, 2010).

The Cenozoic is also a period of significant and relatively rapid climatic and geological changes (Summerhayes, 2015). The biggest geological changes occurred in Eurasia, resulting from its collision with the African, Arabian and Indian plates, that caused the closure of the Tethys Sea and the Alpine - Himalayan orogeny (Trifonov et al., 2012). The climate also changed dramatically during Cenozoic; transitioning Earth from a warm and humid greenhouse state to a cool and dry icehouse; but with a complex interplay of large-scale changes and smaller-scale fluctuations, creating a diverse range of climatic conditions on the local scale (Mudelsee et al., 2014).

Geologic and climatic changes are major drivers of biodiversity evolution (Antonelli et al., 2018; Mayhew et al., 2012; Fritz et al., 2016). Numerous biogeographical studies show distribution patterns and local vicariance and radiation events in the Eurasian biodiversity (e.g. Andreyenkova et al., 2021; Ishii et al., 2023). The majority of these studies focused on terrestrial fauna and flora, but we believe that the best suited model organisms for biogeographical analysis are strict freshwater species. Such hololimnic species are confined to a certain water body unless a physical connection with another water body appears. At the same time, the evolution of the earth’s hydrographic network bases on local geological and climatic conditions. The locally ‘locked’ freshwater fauna thus reflects the environmental past of any region on earth in much more detail than most terrestrial animals or plants (e.g. Bănărescu, 1992; Bohlen et al., 2020a).

Despite this potential, there is a paucity of research on the evolutionary history of Eurasian freshwater biodiversity covering the large geographic and temporal scales. Numerous studies deal with certain geologic or climatic changes during short parts of the Cenozoic (e.g. Kahlke 2014; Ballarin & Li 2018; Gu et al. 2022) or cover a comparably small region of Eurasia (e.g. Delić et al. 2020; Deng et al. 2020; Saito et al. 2018), but there is no comprehensive, overarching analysis that covers the whole of Eurasia, most of the Cenozoic and considers climatic as well as geologic events as driving forces for the evolution of biodiversity.

For the present study we chose the species-rich (> 700 species) freshwater fish family Nemacheilidae (stone loaches) as model; which is omnipresent in almost all water bodies of Europe and Asia (Kottelat, 2012; Dvořák et al., 2022a) and stretches from Spain to Japan in west-east direction and from Russia to indonesia in north-south direction. Most species are small (< 12 cm SL), benthic and inhabit streams and rivers, although some are found in swamps, caves and lakes. Few even occupy extreme environments: Nemacheilidae include the deepest freshwater cave fish in the world (Triplophysa gejiuensis, 400 m below surface) and the fish living at the highest altitude (Triplophysa stolickai, 5200 m a.s.l.) (Chu & Chen 1979, Kottelat 2012).

The aim of the study is to reconstruct the evolutionary history of the family Nemacheilidae as example for the evolution of freshwater fauna across Eurasia and since early Cenozoic. The research concept includes a phylogenetic reconstruction based on molecular genetic data, followed by the dating of cladogenetic events and a biogeographic analysis in order to reconstruct ancestral distribution ranges. The results are compared with the known geological and climatic history of Eurasia to identify the causative factors and processes that have influenced the diversification of Eurasian freshwater fishes.

Results

Major clades within Nemacheilidae

To understand the phylogenetic diversity of the family Nemacheilidae and to identify the largest units within, we analysed sequence data for one mitochondrial (cytochrome b) and up to five nuclear (EGR3, IRBP2, Myosin6, RAG1, Rhodopsin1) genes of 471 specimens representing 279 species and 37 genera from across Eurasia (Fig. 1; Supplementary Material Table S1). This is by far the most comprehensive genetic dataset of Nemacheilidae ever published. The outgroup included 21 species from 9 related families (Fig. 2). We applied maximum likelihood analysis (ML) in IQ-TREE and Bayesian inference (BI) in mrBayes on the concatenated dataset. Besides, the individual unrooted ML gene trees were reconstructed in IQ-TREE and subsequently used as input files for ASTRAL III (Zhang et al., 2018) to infer species tree. All analyses revealed Nemacheilidae as monophyletic group and recovered the presence of six monophyletic major clades within the dataset.

Fig. 1.

Fig. 2.

Each major clade exhibited a distinct distribution (Fig. 3). Overlap of ranges occur in the Burmese region and in Indochina.

Fig. 3.

The most widely distributed is the ‘Northern Clade’ (Fig. 3A), which inhabits the northern parts of Asia and most of Europe from Japan to Spain and southwards to Yunnan. The ‘Southern Clade’ (Fig. 3A) has the second largest distribution and inhabits the Indian subcontinent, the Burmese region and West Asia from the Salween River to southeast Europe; isolated species are found also in Indochina and East Africa (Ethiopia). Another clade is the ‘Burmese Clade’ (Fig. 3B) with a continuous distribution from western Thailand to Northeast India, plus isolated species in the Western Ghats of India and in eastern Thailand and Cambodia. The ‘Sundaic Clade’ (Fig. 3B) includes only the species of the genus Nemacheilus and is distributed across Sundaland (the Great Sunda Islands Borneo, Sumatra and Java and the Malay Peninsula) and in the Indochinese rivers draining to the Gulf of Thailand (Mekong, Mae Khlong and Chao Phraya). The ‘Eastern Clade’ (Fig. 3C) is found along the eastern margin of Asia from the Russian Far East to eastern Borneo. Its disjunct distribution consists of a single species in northeast Borneo, the majority of genera and species in northern Vietnam and southeast China and a single genus on the Korean Peninsula, the Japanese Archipelago, northeast China and the Russian Far East.

The sixth clade (‘Indochinese Clade’, Fig. 3C) has a continuous distribution throughout Thailand, Laos, Cambodia, Vietnam, central Myanmar and southern China, as well as isolated species in the upper/middle Yangtze and on Tioman Island (Malaysia). The distribution areas of the six major clades are not exclusive, but show little overlap in the cases of Northern, Indochinese and Southern Clades, while the Sundaic Clade shares the northern half of its distribution area with the Indochinese Clade and the Eastern Clade even shares its northern part of distribution with the Northern Clade, the middle part with the Indochinese Clade and the southernmost point of occurrence with the Sundaic Clade.

Analyses of the concatenated dataset as well as ASTRAL analysis inferring a species tree from single gene trees revealed a basal dichotomy in the family Nemacheilidae (Fig. 1,2,4). The one branch consists of the Eastern Clade as sister to the pair Indochinese plus Northern Clade. The second branch is made of the Sundaic Clade as sister to the pair Burmese and Southern Clade. All basal nodes and most of tip nodes are supported by high statistical supports in all analyses.

Fig. 4.

Divergence times

The ages of the genealogic events were determined using three fossil records and one geological event as calibration points. The results indicate that Nemacheilidae shared the last common ancestor with other families of Cobitoidea about 48.5 mya (Fig. 4). The first notable split that separated the family into two big branches took place around 39 mya (for details see Fig. 5).

Fig. 5.

Biogeographic reconstructions

The biogeographic reconstruction of Nemacheilidae considered 11 geographic regions within Eurasia and northeast Africa. The results show that Nemacheilidae originated in the region of nowadays Indochina and expanded from there in multiple waves. More details are given in Fig. 6.

Fig. 6.

Discussion

The evolutionary history of Nemacheilidae (Fig. 7)

Early Eocene (56 – 48 mya) - the origin of Nemacheilidae (Fig. 8A)

The family Nemacheilidae shared the last common ancestor with the other Cypriniformes about 48.5 mya (Fig. 2), and our biogeographic analysis identified the region of nowadays Indochina as most likely place of origin of Nemacheilidae (Fig. 6). During the early Eocene, Southeast Asia experienced a warm and perhumid climate, characterised by both highlands and lowlands with numerous rivers (Morley 2018), which must have been ideal conditions for stone loaches.

Fig. 7.

Fig. 8.

Middle to late Eocene (48 – 34 mya) – first range expansions and formation of the first major clades (Fig. 8B,C)

The first detectable cladogenetic event within the family occurred about 39 mya and separated a lineage formed by the Eastern, Northern and Indochinese Clades (Branch A in Fig. 5) from a lineage including the Sundaic, Burmese and Southern Clades (Branch B). Our biogeographic analysis of the dataset revealed Indochina as most parsimony area of the first cladogenetic event; and consequently as region of origin for branch A as well as branch B.

The Eastern Clade was the first of the six present major clades to form: it separated about 35.8 mya from the common ancestor of the Indochinese and Northern Clades, either in the region of Indochina or Vietnam-south China. Taking into consideration its present-day distribution (Fig. 3) we assume that the young Eastern Clade spread mainly along the eastern coastline of Asia southwards to Borneo and northwards to northeast China. Geologic data support this hypothesis: During the Eocene, Asia was much smaller than today and had a rather flat relief with old mountains mainly along its eastern coast due to the collision with the Pacific Plate since 190 mya (Seton et al., 2012; Wang, 2004), while the upfolding of east-west mountain ridges had just started. The generally warm and humid climate during Eocene (Morley 2018; Westerhold et al. 2020) promoted the range expansion of the Eastern Clade, while the ancient eastern mountain ranges directed the range expansion. It can be assumed that shortly after the range expansion the distribution of the Eastern Clade was more or less continuous and not separated into several disjunct areas as nowadays.

The ancestral Indochinese + Northern Clades also spread northwards, but their present-day distribution (Fig. 3) suggests that they took a route through Central Indochina and China. Most likely, this first geographic expansion led the ancestral Indochinese+Northern Clades northward to at least the Tibetan Protoplateau, which during Eocene epoch had its northern boundary in the region of the Tanggula Shan (Wang et al., 2008). The presence of the two protoclades on the early Tibetan Plateau is indicated by the presence of some of their most basal members in the Yangtze basin and Yunnan: the genus Homatula for the Indochinese Clade and several species of Triplophysa for the Northern Clade (Figs. 1,5).

Oligocene (34 – 23 mya) – range expansions, extinctions and the formation of the remaining major clades (Fig. 8D,E)

The Eocene - Oligocene transition brought a period of significant global cooling (Hutchinson et al 2021) as well as a change in the precipitation regime of Central Asia (Clift & Webb 2018; Ye et al. 2022) due to the ongoing uplift of the Himalayan Mountains (Ding et al. 2017). A substantial decrease of rainfall in Central Asia caused an aridification from 34 – 23 mya (Bosboom et al. 2014) extending across nearly all of Oligocene Asia from the western to the eastern coast (Wang, 2004). The distribution area of the Eastern Clade became fragmented during this period (Fig. 3).

The Central Asian arid zone also split the Indochinese – Northern protoclade into two part that evolved into two separate major clades, the Indochinese Clade and the Northern Clade. Moreover, the distribution area of the Indochinese Clade was split and separated the southernmost branch (Speonectes and Tuberoschistura) from the northern branch. The latter presently inhabits all of Indochina and southern China including the middle Yangtze basin (Fig. 3), indicating that the Oligocene refuge of the northern branch of the Indochinese Clade was located southeast of the Tibetan Plateau. In contrast, the early Northern Clade found its refuge further north on the Tibetan Protoplateau, the only region that can explain the further history and the present distribution of the clade. By late Oligocene, the climate in Central Asia returned to a warmer and moister state (Zachos et al. 2001) and the eastern Kunlun Shan Mountains rapidly thrust upward as an expansion of the northern Tibetan Plateau (Li et al. 2021). These favourable conditions coincided with a major split within the Northern Clade at 20 mya (Fig. 5), when the proto-Barbatula clade colonised the Kunlun Mountains, while the proto Triplophysa branch found its main distribution rather east of the Tibetan Plateau.

Different geologic events drove the evolution of branch B: the common ancestor of the Burmese and Southern Clades split at 32.5 mya in Indochina from the Sundaic Clade. In contrast to the general trend among Nemacheilidae, the species of the Sundaic Clade (genus Nemacheilus) prefer sandy bottom over gravel bottom, making lowlands regions of Sundaland and Indochina their ideal habitat. There are no indications that the Sundaic Clade ever left the region of present-day Southeast Asia.

The separation of the ancestral Burmese and Southern Clades from the Sundaic Clade occurred because the former left Southeast Asia by colonising the West Burma Terrane, a microplate that was pushed northward by the Indian plate (Morley et al. 2021) and made temporary contact with Sundaland around Eocene-early Oligocene (Westerweel et al. 2019). Around 30 mya, it occasionally made contact with northeast India (Morley et al. 2021), allowing fishes from the West Burma Terrane to colonise the Indian subcontinent. During the Oligocene, the Indian plate and the West Burma Terrane both collided with mainland Asia. The former presently forms the Indian subcontinent, the later Western Burma (Westerweel et al. 2020). Their parallel collision with mainland Asia caused the uplift of the north-south directed Indo-Burmese mountain range (Morley et al. 2020), providing a physical barrier between the Brahmaputra basin on the Indian side from the Chindwin-Irrawaddy basin on the Burmese side, promoting the separation of the two faunas into the Burmese Clade and the Southern Clade.

Miocene (23 – 5.3 mya) – rapid diversification (Fig. 8F,G)

The Miocene is generally an epoch of cooling, with temperatures decreasing about 5.5° C to conditions similar to present day (Zachos et al. 2001). But temperatures also underwent strong fluctuations, leading to exceptionally warm periods like the Middle Miocene Climatic Optimum (17 – 14.5 mya) as well as to drastic cooling events like the Middle Miocene Climatic transition (14 mya) (Morley 2018). Additionally, the Asian monsoon regime strengthened from the early Miocene (Clift & Webb 2018), increasing precipitation and altering the distribution of climatic regions across Asia (Sun & Wang 2005). Moreover, during the Miocene the Himalayan Mountains experienced their main uplift (Ding et al. 2017) and the upfolding area of Central Asia mountain ranges expanded northwards, giving rise to the Tian Shan and Altai mountains (Nissen et al 2009; Liu et al 2023).

During this complex geological and climatic scenario, the fish family Nemacheilidae underwent a massive diversification. In our phylogenetic hypothesis the family comprised 16 lineages at the beginning of Miocene and 177 lineages by the end of the epoch. The radiation occurred across all six major clades, albeit to a lesser extent in the Eastern Clade. Interestingly, it manifested similarly in the region of the Indian monsoon and the East Asian monsoon. In both regions, major radiations occurred during the Middle Miocene Climatic Optimum, including the radiation of the Barbatula lineage in the Northern Clade, of the Indochinese Clade, of all Nemacheilus except N. binotatus and N. ornatus and of six branches in the Southern Clade. It appears that under increased temperature and precipitation freshwater fishes flourished, across all of Asia and under different monsoon systems. One of the radiations that began during the Middle Miocene Climate Optimum involves Afronemacheilus, the sole nemacheilid genus found in Africa (Prokofiev & Golubtsov 2013). Our results confirm that the ancestral Afronemacheilus separated from the Southern Clade about 18 - 17 mya. It could have utilised two potential pathways to Africa: one across the Arabic Peninsula and via the southern Danakil Isthmus, a broad landbridge that stretched from Djibouti to Aden from late Miocene to Pliocene (Bosworth et al. 2005; Steward & Murray 2017). The second pathway would lead through the Levant and via Sinai into the Nile (‘Gomphotherium landbridge’ Harzhauser et al 2007), and subsequently upstream into the Ethiopian highlands. Unfortunately, the nearly complete extinction of freshwater fishes from the Arabic Peninsula and from the Nile during Pleistocene aridification (Kotwicki & Salaimani 2009; Leplongeon 2021) complicates the reconstruction of its pathway to Africa.

In the Northern Clade, the Triplophysa lineage flourished during the Miocene. Its present wide distribution in southern and Central China suggests that it expanded from its Oligocene refuge on the eastern Part of the Tibetan Plateau in an eastern direction (Fig. 3). In contrast, the Barbatula branch expanded from its refuge on the western Tibetan Plateau northwards, traversing the uplifting Tian Shan and Altai mountains. This expansion can be traced from the geographic position of the most basal species (T. labiata-Tian Shan, B. altayensis – Altai), and then eastwards (B. compressirostris + D. dzebuadzei – west Mongolia, B. cf. toni – east Mongolia; tip species in Korea (B. indet Korea) and Japan (B. oreas).

During late Miocene, Barbatula also expanded westwards through Siberia and colonised Europe between 5.6 mya (separation of the European taxa from the Mongolian sister taxon) and 3.0 mya (oldest split within European Barbatula). During this cool and dry period, first glaciations occurred on the northern hemisphere (Zachos et al. 2001), and in Siberia the mean temperature of the coldest month dropped below 0° for the first time (Popova et al. 2012). Freezing of the northern margin of Eurasia impeded the Siberian rivers (which flow from south to north) from discharging their water into the Northern Ocean and produce extensive wetlands on the shallow relief of Siberia. These wetlands served as pathways for freshwater fishes from East Asia to European river systems during Pleistocene (Kalous et al. 2012), and we assume that the same mechanism brought Barbatula to Europe during Pliocene.

Pliocene and Pleistocene (5.3 – 0.01 mya) - Ice ages

During Pliocene and Pleistocene global cooling continued until massive glaciations repeatedly covered the north of Eurasia from early Pliocene onwards (de Schepper et al. 2014). However, the Pliocene also brought a much stronger seasonality to Eurasia, especially in the form of cold winters (Shen et al. 2018). Despite these significant changes in climatic conditions, the evolution of Nemacheilidae seems not to have experienced any dramatic negative impact, but reveals mainly local speciation events (e.g. in the genera Traccatichthys and Barbatula). We know from more detailed studies on Nemacheilidae (Bohlen et al 2020a,b; Šlechtová et al. 2021) that during Pliocene and Pleistocene also geographic shifts as well as expansions and disappearances occurred, but since these happened on rather local scale they were too insignificant to reflect in the present intercontinental dataset.

Conclusions

• - We reconstruct the evolutionary history of the Eurasian fauna during most of Cenozoic era on the example of the freshwater fish family Nemacheilidae. The extensive dataset (471 specimens, 279 species, 6 genes) from across Eurasia and the low dispersal capacity of the model resulted in high resolution of the phylogenetical and biogeographic analysis.

• - Molecular phylogeny uncovered six major clades within the family. Dating of cladogenetic events and ancestral range estimation traced the origin of Nemacheilidae to Indochina around 48 million years ago.

• - Major expansions were facilitated by tectonic connections, favourable climatic conditions, and orogenic processes. Conversely, aridification emerged as the primary cause of extinction events.

• - Our study represents the first comprehensive and holistic reconstruction of the evolution of Eurasian freshwater biodiversity on a continental scale and through deep geological time.

Methods

Sampling and Laboratory procedures

471 samples of more than 250 species from 37 described genera of the family Nemacheilidae were examined, including the set of sequences of 364 specimens analysed in our laboratory and sequences of 107 specimens obtained from GenBank. Specimens used for the present study were fixed and stored in 96% ethanol. For more details about samples, sampling sites and GenBank accession numbers see Supplementary material Table 1.

Total genomic DNA was isolated using the DNeasy Blood & Tissue kit (Qiagen) following the manufacturer’s instructions. One mitochondrial gene (cytochrome b) and five nuclear genes (RAG1, IRBP2, MYH6, RH1 and EGR3) were selected. For list of primers see Supplementary material Table 2, details about PCRs settings are described in Dvořák et al., 2022; Chen et al., 2003 and Chen et al., 2008.

Sequencing of the PCR products was performed in Laboratory of Fish Genetics on ABI Prism 3130 GA or via sequencing service in Macrogen Europe BV (Amsterdam, Netherlands).

Phylogenetic analyses

Chromatograms were checked and assembled in the SeqMan II module of the DNA Star software package (LASERGENE). Single gene alignments were done in BioEdit (Hall, 1999) with use of ClustalW (Larkin et al., 2007) multiple alignment algorithm. TheThe datasets were concatenated in PhyloSuite v1.2.2 (Zhang et al., 2020).

The best-fit substitution models and partitioning schemes (Supplementary material Table 3) were estimated using Partition Finder 2 (Lanfear et al., 2016) implemented in PhyloSuite 1.2.2 based on the corrected Akaike Information Criterion (AICc).

The 471 specimen phylogenetic trees were inferred from six loci concatenated dataset using the maximum likelihood (ML) and Bayesian inference (BI). The ML analyses were performed using IQ-TREE (Nguyen et al., 2015) implemented in PhyloSuite. The best-fit evolutionary model for each codon partition was automatically determined by ModelFinder (Kalyaanamoorthy, et al. 2017). The node support values were obtained with 5000 ultrafast bootstrap replicates (UFBoot) (Hoang et al., 2018). For the BI analyses we used MrBayes 3.2.7a (Ronquist & Huelsenbeck, 2003) on CIPRES Science Gateway (Miller et al., 2010). The datasets were partitioned into genes and codon positions and analyses were performed in two parallel runs of 10-20 million generations with 8 Metropolis Coupled Markov Chains Monte Carlo (MCMCMC). The relative burnin of 25% was used and from the remaining trees a 50% majority rule consensus trees were constructed.

For reconstructing of the species tree we used the Accurate Species TRee ALgorithm (ASTRAL III) (Zhang et al., 2018). For ASTRAL we used unrooted single gene ML trees reconstructed in IQ-TREE implemented in PhyloSuite. We have used IQ-TREE 2 (Minh at al., 2020) to calculate gene concordance factor (gCF, indicating the percentage of decisive gene trees containing particular branch) and site concordance factor (sCF, indicating the percentage of decisive alignment sites that support a branch in the reference tree) for every branch in the ASTRAL species tree. Distribution maps in Fig. 2 were designed using HYDROSHEDS (Lehner & Grill 2013).

Divergence time estimations and ancestral range reconstruction

The ages of cladogenetic events were estimated in BEAST 2.6.4 (Bouckaert et al., 2014) via CIPRES Science Gateway with use of four calibration points. The first calibration point is based on the oldest known fossil of the family Catostomidae, Wilsonium brevipinne, from early Eocene (Liu, 2021), approximately 56-48 my old. The second calibration point is derived from the only known nemacheilid fossil record, Triplophysa opinata from Kyrgyzstan from middle-upper Miocene (16.0 to 5.3 mya) (Prokofiev, 2007; Böhme and Ilg, 2003). Third calibration point is based on the oldest known fossil of the genus Cobitis, C. naningensis, from early to middle Oligocene in southern China, 34-28 my old (Chen et al., 2015). The fourth calibration point, based on the river history of the southern Korean Peninsula, dates the separation between the species Cobitis tetralineata and C. lutheri to 2.5-3.5 my (Kwan et al. 2014). For the fossil-based calibration points we have used a fossil calibration prior implemented in CladeAge (Matschiner et al., 2017), the fourth calibration point was set to a uniform prior (2.5 – 3.5). The calibration points are depicted in Figs. 2 and 5.

For the final analyses, the partitions were unlinked and assigned the estimated evolutionary models. We used relaxed lognormal molecular clock and birth-death prior. Given the large and complex nature of our dataset, and in light of several preliminary analyses where MCMC chains did not provide satisfactory ESSs even after 1 billion iterations, made the decision to streamline the dataset by retaining only one specimen per species. To enhance the performance of the analysis, we also implemented parallel tempering using the CoupledMCMC package (Müller & Bouckaert, 2019, 2020). The analysis was configured with four parallel chains of 6×10⁸ generations, with resampling every 1000. Tree and parameter sampling intervals were set to 60.000. The effective sampling sizes (ESS) for all parameters were assessed in Tracer 1.7.1 (Rambaut et al., 2018). A maximum clade credibility (MCC) tree was built in TreeAnnotator 2.6.0 (Rambaut & Drummond, 2010) after discarding the first 25% of trees. The final trees were visualised in FigTree 1.4.4 (Rambaut, 2019). ML and BI were performed with both full as well as the reduced datasets.

The biogeographical analysis was conducted using the BioGeoBEARS package (Matzke, 2013) implemented in RASP 4.0 (Reconstruct Ancestral Stage in Phylogenies, Yu et al., 2015). We set 11 biogeographic regions according to drainage areas separated by major mountain ridges (Ural, Caucasus. Himalayas, Tenasserim Ridge, Indo-Burmese Range, Annamite Cordilliera, Balkan, Hindukush) or saltwaser straits (Red Sea) (see Fig. 4A). For the analysis, we utilized the final MCC tree and 100 post burnin trees from the previous BEAST analysis. A maximum of 2-unit areas was allowed to be present in every reconstructed ancestral range. The comparison of six biogeographic models (DEC, DEC+J, DIVALIKE, DIVALIKE+J, BAYAREALIKE, BAYAREALIKE+J) suggested that the DEC+J is the best-suited model for the given dataset (as indicated by AICc weight value, Supplementary material Table 4). Likelihood Ratio Test (LRT) p-values, comparing models with and without the J factor, indicated that the addition of the J factor significantly influences the model likelihood. Considering recent criticisms of DEC and DEC+J (Ree & Sanmartin, 2018), we opted to conduct the analysis not only with DEC+J but also with other +J models. Comparison of the results from all applied models revealed only minimal differences.

Acknowledgements

VŠ, TD, VŠ and JB would like to thank R. Hoyer, J. Kühne, F.F. Kullander (NRM), Jan Kuszniersz, K. Lim (ZRC), H. Linke, A. Nolte, G. Ott, M. Reichard, I. Seidel, S. Somadee, H.H. Tan (ZRC), T.K. Toe, K. Udomritthiruji, T. Win and U. Wolf for their help with obtaining specimens for this study and Š. Pelikanová for help with laboratory work. VŠ and JB were funded by grants 206/08/0637 and 19-18453S of the Grant Agency of the Czech Republic. BL and AG greatly acknowledge S. Cherenkov for his help in sampling. BL was supported by Russian Science Foundation, grant no. 24-44-20019.

Supplementary information

Table S1.

Table S2.

Table S3.

Table S4.