Phylogenetic analyses and species delimitation
We analysed 104 individuals from the studied complex collected at 64 localities (Fig. 1). In order to explore the phylogenetic position of the focal taxa, we compiled a dataset with additional 37 species (details in Supplementary Information, Table S1) covering all main Niphargus lineages22, 23. Bayesian and Maximum Likelihood analyses yielded identical tree topologies (Fig. 2, Fig. S1). The studied group of taxa was not monophyletic but split into four major non-sister lineages. The first three lineages represent the nominal N. stygius stygius from West Slovenia, the nominal N. s. karamani from North-Eastern Slovenia and the nominal N. s. kenki distributed between Central and Eastern Slovenia along the Sava River. All remaining species are nested within a strongly supported clade together with already described and morphologically distinct species N. illidzensis, N. slovenicus, N. dalmatinus, N. vinodolensis. N. zagrebensis and N. elegans (Fig. 2, Fig. S1). This last clade is distributed between the lowlands of northern Italy, the northern Dinaric Karst and the western margins of the Pannonian Plain.
Species hypotheses were inferred using uni- and multi-locus tree-based24, 25 species delimitation methods, and unilocus distance-based26 methods. We first ran the unilocus Poisson Tree Processes (PTP)24 delimitation procedure. The analyses were performed on a COI alignment counting 117 sequences belonging to herein studied and 37 other species (Supplementary Information 1, Table S1). The results suggest that the studied species complex counts 16 distinct species. Nominal N. s. brachytelson, N. s. novomestanus, N. s. kenki and N. s. karamani each split into two distinct species, nominal N. s. podpecanus into three, and N. s. likanus into four species. Only nominal N. s. stygius remained a well-supported single species. In addition to the focal species complex, the analyses unveiled another nominal taxon hiding additional cryptic diversity, i.e., the nominal N. zagrebensis turned out to consist of two morphologically similar species. Next, we applied a multilocus coalescence delimitation method using Bayesian Phylogenetics & Phylogeography 3.1 (BPP)25. As the N. stygius s. lat. species complex turned out to be non-monophyletic, we ran this analysis separately for each of three lineages that tentatively included additional species. Multilocus species delimitation confirmed the results of the unilocus analysis, regardless of the population size and root age priors used (Fig. 2, Table S3).
In order to assure robustness of species hypotheses, we supplemented the tree-based species delimitations by a distance-based approach. We tested whether COI sequences of the putative species diverged more than 16% in patristic and 4% in Kimura-Two-Parametric (K2P) distances. The two thresholds are conservative estimates of lower species boundaries, determined empirically. The 16% patristic distance corresponds to the divergence among 276 well-defined crustacean species26, whereas a 4% K2P distance indicates the divergence beyond which interbreeding among Gammarus amphipods becomes unlikely27. Pairwise comparisons of the delimited species indicate high molecular divergence between most species pairs. All species pairs exceeded the 4% K2P threshold, and all but four species pairs exceeded the 16% patristic distance threshold (Table S4, Fig. 3). However, even in these four species pairs, interspecific distances substantially exceed intraspecific ones (Fig. 3).
As a final criterion, we used sympatric occurrence of cryptic phyletic lineages as a pointer to reproductive isolation and hence distinct species even under the biological species concept (Figs 1 and 2). We found that out of the 16 delimited species ten pairs had partially overlapping ranges (Fig. 1). In five caves, we encountered species pairs in syntopy (i.e., three species pairs, of these one pair syntopically co-occurred three times). These observations we treated as natural crossing experiments, and they revealed no sign of gene flow or introgression between species. Among one mitochondrial and three nuclear loci all haplotype were consistently species-specific.
Four boundary cases, which diverged less than 16% in their patristic COI distances, are allopatric. However, a close examination unveils that in three of these species pairs divergences arose over short geographic distances (10–20 km, Fig. 3). These species pairs were found within the same catchments, indicating that the observed divergences and uniqueness of lineages may be explained by biological reproductive barriers rather than limited dispersal. Less clear is the case of the least diverged species pair (N. kenki complex; 9% patristic and 5.5% of K2P distance) separated by 80 km. Although all delimitation methods but one suggest this last pair of populations should be treated as separate species, the data are insufficient to rule out alternative explanations, such as isolation by distance. For this reason we decided to postpone taxonomic decisions about distinct N. kenki populations until additional data clarify their status.
Impact of cryptic species on species richness
The results of the species delimitation procedures substantially affected species richness measures. The species number has risen by a factor of two and up to a factor of 15, depending on whether traditional subspecies are treated as species28 or not21, 29. We analysed the effects of change in species richness at two geographic scales that are relevant for conservation: at single sites and nation-wide. As the complex is spread across the Danube and Adriatic basins, we investigated whether species richness equally increased in the two basins.
At the level of individual caves, the effect of increased species richness was negligible. Species richness has increased by one species in five caves where syntopic occurrences of two species were detected. On a nation-wide scale, the number of species increased by a factor of 12 in Slovenia and by a factor of five in Croatia. Interestingly, if subspecies were treated as phylogenetic species28, the increase in species richness in both countries would be approximately two fold (1.9 for Slovenia and 1.7 for Croatia). Increase in species richness in the two basins was strongly asymmetric. All newly discovered cryptic species were detected in the Danube basin.
Changes in range size – a proxy of species vulnerability
Changes in range sizes were remarkable. The range of the nominal N. stygius sensu S. Karaman, including all its subspecies, measured over 20,000 km2, spanning East Italy, Slovenia and West Croatia. This range is unusually large for an aquatic species living in subterranean environment12. It is further unusual in that it covers the Adriatic and the Black Sea Drainage, as well as four major European biogeographic regions: the Alpine, the Mediterranean, the Dinaric and the Pannonian30. By contrast, ranges of cryptic species within the complex are much narrower. Five species are known from a single locality only and two species are known from two localities. Among species known from more than three localities, three have ranges smaller than 500 km2, and six have ranges between 500 and 5000 km2. As for comparison, according to the IUCN Red List criteria, range sizes below 5000 km2 indicate species are endangered31.
Evolutionary distinctness of cryptic species
Another measure relevant for conservation is evolutionary distinctness (ED), a metrics that estimates evolutionary uniqueness of a species as compared to its congeners32, 33. The measure may be related also to a species’ function in the ecosystem (see Discussion)34. ED is a measure of a species’ terminal branch-length, corrected for the species richness of the containing clade; i.e., species on longer branches with fewer congeners receives higher ED value than species within recent, species rich radiations32. It is measured either in millions of years or in the amount of nucleotide substitutions accumulated along the branch of a phylogenetic tree.
The evolutionary distinctness of N. stygius as traditionally conceived was 0.087 nucleotide substitutions per site. As cryptic species belonged to four independent clades within the wider phylogeny of Niphargus, their ED is not diminished by splitting. The ED values of the herein studied species of the N. stygius complex differed from each other by a factor of 2.4 (0.037–0.087, Table 1). As these values are little informative perse, we compared them to the values calculated for the other, non-cryptic species included in our phylogenetic analysis. About one third of the species attained relatively high (>0.075), and another third low (<0.045) ED values (Fig. 4). Overall, the ED values of the studied complex were only slightly, but not significantly lower than in the non-cryptic species from the rest of the phylogeny (Mann-Whitney U test, P = 0.092).
Species diagnoses and names
We showed that, according to different criteria, the 15 delimited genetic lineages represent fully evolved species that likely meet even the most rigorous criteria of biological species35. Moreover, because of their small ranges, exceptionally high endemism, and their evolutionary distinctness, they are highly relevant for conservation. We therefore undertook the necessary formal steps to name and diagnose them as species under the provisions of The International Code of Zoological Nomenclature (The Code). We raise seven existing subspecies to species rank, and diagnose nine species new to science.
The morphological diagnosis for the entire complex as provided by S. Karaman21 remains valid and is available in Supplementary Information 3. For diagnosing the species, we followed recent recommendations17, 18 (details in Material and Methods); the analyses were made in CAOS36. All species were diagnosed using unique combination of character attributes at four genes. The number of diagnostic characters varied between 0 and 156 per marker (Table 2). Specimens from localities from which subspecies were described were considered as representatives of already named taxa. An intriguing problem was to assign the name N. podpecanus to the right species. In its type locality, the cave Podpeška jama, we found two co-occurring, genetically distinct but morphologically indistinguishable species. The type series is approximately 70 years old and no longer contains any useful DNA. However, in his original work, S. Karaman21 reported that the species was found also in the vicinity of the city of Kočevje. In order to maximize compatibility between the bibliographic record and our new taxonomic conclusions, we assigned the name N. podpecanus to the one species that we found also in Kočevje. Following Article 75 of the Code, we erected a neotype for N. podpecanus (Table 2).
Species names, etymology, voucher numbers and information about type repositories are available in Table 2 and in Supplementary Information 3. The specimens and samples are available for further exploration; alignments and used for diagnosing the species are deposited at Dryad (doi:10.5061/dryad.6rs6q). In order to make descriptions compliant with The International Code of Zoological Nomenclature, the species were registered in Zoobank and World Register of Marine Species (Table 2, Supplementary Information 3).
Not every article in a journal is considered primary research and therefore "citable", this chart shows the ratio of a journal's articles including substantial research (research articles, conference papers and reviews) in three year windows vs. those documents other than research articles, reviews and conference papers.