Aim: We tested the hypothesis that historical marine regression-transgression cycles shaped the distribution patterns of subterranean amphipods through repeated cycles of dispersal and vicariance against the hypothesis that subterranean amphipods colonized both sides of the Adriatic Sea independently.

Location: Western Balkan Peninsula, Adriatic Sea Islands and Apennine Peninsula, Europe

Taxon: A clade of freshwater subterranean amphipods, genus Niphargus (Crustacea: Amphipoda).

Methods: The taxonomic structure of the studied clade was revised using unilocus species delimitations. Timeframe of cladogenetic events was inferred using multilocus time-calibrated phylogeny and compared to the main regression-transgression events in Miocene and Pleistocene. The geographic origin of the studied clade, species’ range expansions and contractions, as well as vicariant events were assessed through modelling of historical biogeography.

Results: Subterranean amphipods of the genus Niphargus, found on both sides of the Adriatic Sea, form a monophyletic clade. The reconstructions of ancestral ranges suggest that the clade emerged in the Balkan Peninsula, three times independently dispersed to the Apennine Peninsula and once back to the Balkans. Adriatic Islands were colonized multiple times, predominantly from the Balkan Peninsula. The dispersal-vicariance events correspond to historical regression-transgression cycles in Miocene and Pleistocene.

Main conclusions: Marine regression-transgression cycles shaped distribution patterns of subterranean amphipods, while the alternative hypothesis received no support. The subterranean faunas apparently well reflect older biogeographic events.

Taxon sampling and DNA isolation

Taxon sampling aimed to include a wide range of Niphargus species, distributed from Ireland to Iran, and to comprehensively sample species distributed on the both sides of the Adriatic Sea. We sampled altogether 214 localities and totally 494 individuals, and stored them in 96% ethanol. 

Each specimen was identified to the closest morphologically resembling morphospecies using original species descriptions. 

Taxonomic structure and molecular species delimitation

In order to account for incomplete taxonomy and the possible presence of morphologically cryptic species, we revised the taxonomic structure of the study dataset prior to the main analyses. We assembled a dataset containing altogether 463 COI sequences of Niphargus and genera nested within it (Carinurella, Haploginglymus, Niphargobates, Pontoniphargus) The sequences were acquired in this study or retrieved from GenBank. Putative species structure, represented by the so called molecular taxonomic operational units (hereafter MOTUs) was inferred using distance- and tree-based delimitation methods, automatic bridge gap discovery (hereafter ABGD) (Puillandre, Lambert, Brouillet, & Achaz, 2012) and Poisson tree processes (hereafter PTP) (Zhang, Kapli, Pavlidis, & Stamatakis, 2013).

Phylogenetic analyses

Phylogenetic structure of the studied dataset was inferred using 195 Niphargidae species (Niphargus and the genera nested within) rooted with five species of Pseudoniphargus, which was shown to be a sister group in two previous studies (Jurado-Rivera et al., 2017; Copilaş-Ciocianu et al., 2020). The dataset included 169 molecular operational taxonomic units (MOTUs) delimited in the PTP analysis (see above and Results section) supplemented with 31 described species. Phylogenetic relationships were reconstructed using two alternative approaches: i) maximum likelihood with partition-specific setting and ultrafast bootstrapping in IQ-tree 1.6.7 (Nguyen, Schmidt, Von Haeseler, & Minh, 2015), and ii) Bayesian inference with partition-specific settings in MrBayes 3.2.6 (Ronquist et al., 2012).

Estimation of divergence times

In order to define the timeframe of splits between the lineages distributed in the Dinaric Karst and the Apennine Peninsula, we reconstructed a time-calibrated multilocus phylogeny using the package BEAST 2.5.1 (Bouckaert et al. 2018). The dataset comprising 200 individuals and four molecular markers was analysed using partitioning schemes with distinct parameters settings, following best-fit models of evolution proposed by bModelTest 1.2.1 (Bouckaert & Drummond, 2017). Substitution and clock models were unlinked for all partitions, while the tree partitions were linked to generate a single dated phylogeny. Birth death and Yule speciation tree priors and their marginal likelihoods were compared in BEAST’s extension Path Sampler. Based on marginal likelihood, we used a Yule speciation tree prior for the final analyses (results not shown). In order to account for lineage-specific rate heterogeneity, we employed a lognormal relaxed clock (Drummond, Ho, Phillips, & Rambaut, 2006). 

Modelling past dispersal, vicariance and reconstruction of ancestral ranges

In order to elucidate the historical biogeography of the main Transadriatic clade, we modelled dispersal, extinction, vicariance, cladogenesis and ancestral ranges (Fig. 1) using package BioGeoBEARS 1.1.2 (Matzke, 2013) in R 3.5.3 (R Core Team, 2019).