A hotspot of groundwater amphipod diversity on a crossroad of evolutionary radiations

Groundwater harbours an exceptional fauna and provides invaluable ecosystem services, yet is among the least explored and consequently least protected ecosystems. Successful protection of its biodiversity depends on complete species inventories, knowledge of species spatial distribution, and quantification of biodiversity patterns, as well as disentanglement of the processes that shaped biodiversity patterns. We studied the hyper‐speciose amphipod genus Niphargus as a model system within a global subterranean biodiversity hotspot. We linked the biodiversity patterns with possible underlying processes and discuss the needs to include information on different origins of biodiversity into conservation approaches.


| INTRODUC TI ON
Groundwater is the largest reservoir of liquid freshwater globally that is of direct and major interest to humans (Gibert & Deharveng, 2002;Mammola et al., 2019). Despite its vastness, groundwater is also one of the most endangered habitats. Habitat loss, overexploitation, contamination and climate change threaten its inhabitants (Mammola et al., 2019). Groundwater fauna is unique in many aspects. It comprises phylogenetic and geographic relicts (Humphreys, 2000), disproportionally high number of endemic species (Bregović et al., 2019) and plays a key role in ecosystem services (Griebler & Avramov, 2015). However, conservation biologists have largely neglected groundwater biodiversity, even in comparatively well-studied regions such as Europe and North America (Mammola et al., 2019). Species inventories and mapping of groundwater biodiversity have been lagging behind initiatives on surface systems (Ficetola et al., 2019). Among the main reasons for this are limited access to groundwater hindering sampling , deficiency in taxonomic expertise combined with high morphological similarity (i.e., cryptic species) , and limited availability of robust phylogenies. Consequently, and due to contemporary human activities, many species may be doomed before even being discovered (Niemiller et al., 2013). To reverse this trend, increasing scientific knowledge is crucially needed for the protection of groundwater (Wynne et al., 2021).
Successful conservation management must be based on solid scientific evidence, which itself depends on sufficient empirical knowledge and understanding of the distribution of biodiversity.
The first step in assessment of biodiversity are species inventories at an appropriate spatial and taxonomic resolution (Guralnick et al., 2018). The second step is quantification of biodiversity itself. Because of phylogenetic relatedness and ecological specialization, species cannot be treated as independent or functionally equivalent units. For this reason, quantification of biodiversity has been increasingly relying on combining complementary metrics, such as species richness and phylogenetic diversity (Faith, 1992).
Phylogenetic diversity arguably indirectly captures also other aspects of biodiversity, such as functional diversity (Cadotte et al., 2011;Isaac et al., 2007;Mazel et al., 2018) and has been incorporated into identification of hotspots and assessments of processes that contributed to high biodiversity (Devictor et al., 2010;Isaac et al., 2007;Purvis et al., 2000). In groundwater, biodiversity patterns were hitherto mainly inferred using species richness (Bregović et al., 2019;Eme et al., 2018;Zagmajster et al., 2014), while other metrics of biodiversity have not been applied yet. Eventually, the assessment and understanding of patterns of biodiversity using different metric may also reveal the underlying processes that shaped these patterns.
Deficiency in knowledge and incomplete application of biodiversity measures is especially problematic in areas of exceptional biodiversity, as it may lead to incorrect yet highly influential conclusions with respect to conservation. Such an area is the Western Balkans, which covers 166,000 km 2 in south-eastern Europe, including the Dinaric Karst and the Eastern parts of the Southern Calcareous Alps.
The region is a global biodiversity hotspot for surface and subterranean ecosystems (Myers et al., 2000;Sket, 2012;Zagmajster et al., 2014). Importantly, groundwater biodiversity of this region is locally highly endangered due to planned construction of dams, roads and other infrastructure, and agricultural developments, which lead to habitat destruction (Schwarz, 2015). Despite a rich speleobiological history in the region, much of the groundwater diversity is still being discovered (Zagmajster et al., 2014), and its biodiversity has not yet been assessed using metrics beyond simple species counts.
Here, we studied groundwater biodiversity patterns within the Western Balkans using the subterranean amphipods of the genus Niphargus. With more than 400 known species and still a high number of undescribed species, Niphargus is an exceptionally diversified genus that importantly co-shapes biodiversity patterns in aquatic subterranean habitats of the Western Palearctic Fišer, 2019;Horton et al., 2021;Zagmajster et al., 2014).
Phylogenetic relationships within the genus are relatively well understood, and the evolutionary history of the genus precedes the origin of the Western Balkans as a geologic unit .
Extensive field explorations over the past decade resulted in a large number of fresh Niphargus samples. This offered an opportunity to conduct a study entirely based on molecular data to assess biodiversity patterns of this region using species richness and phylogenetic diversity. We used these metrics to link biodiversity patterns with historical processes that shaped the patterns, and to discuss how these processes should be considered in conservation planning.
We first assessed diversity of Niphargus using molecular unilocus delimitations (Eme et al., 2018) and analysed the spatial and taxonomic coverage of the molecular dataset by comparison with nominally described and published species distribution records (Bregović et al., 2019;Zagmajster et al., 2014). Then we constructed a timecalibrated phylogeny to quantify phylogenetic diversity. A comparison of species richness patterns with phylogenetic diversity allowed us to assess whether species in the richest cells assembled through local diversification (low phylogenetic diversity) or colonization of unrelated species (high phylogenetic diversity). Finally, we discuss the importance of different historical processes in conservation biology.

K E Y W O R D S
amphipods, evolutionary radiations, groundwater, Niphargus, phylogenetic diversity, species richness, Western Balkans 2 | ME THODS

| Overview of the dataset
We examined the collection of Niphargus from SubBio Lab, Department for Biology, Biotechnical Faculty, University of Ljubljana, Slovenia, containing >11,000 individuals, and being one of the largest collections of Niphargus worldwide. We identified the specimens using morphological diagnosis and cross-examined them with existent molecular and morphological data in SubBio Database.
We then selected 864 specimens that were sequenced for the standard barcoding marker cytochrome C oxidase subunit I (COI), to cover all existent localities and as much morphological variability as possible. Combined with already published genetic material from the Western Balkans and 280 representatives of species outside the study area, the final COI delimitation dataset counted 1492 specimens. After cross-examination of delimitations and quality checking, the dataset used for diversity patterns assessment was comprised of 1212 individuals that were aligned to 245 MOTUs, present in altogether 598 localities within the study area, of which 297 localities did not have previous published molecular data for Niphargus (see Appendix S1 in Supporting Information). Although this is the most extensive dataset with the highest coverage for any subterranean aquatic group in this region and well beyond (Zagmajster et al., 2014), it is based on a broad number of studies and datasets, as is common for larger faunistic databases. In order to check for consistency in sampling completeness within the studied region, we thus plotted species accumulation curves (Gotelli & Colwell, 2001) and rarefaction and extrapolation sampling curves (Chao et al., 2014) for the two main subregions, namely the north-west and south-east (see Appendix S5, Figure S5.1). While none of the curves reached a plateau, the result suggests relative homogeneity in sampling. We conclude that results and conclusions are thus comparable and conservative, as, if any, the south-eastern part is closer to the saturation of sampled biodiversity.

| Molecular analysis and delimitation procedures
We extracted genomic DNA from 864 specimens using MagMAX DNA Multi-Sample Kit (Thermo Fisher Scientific). Oligonucleotide primers and amplification protocols were the same as in . All nucleotide sequences were obtained by Macrogen Europe laboratory (Amsterdam, Netherlands), using amplification primers and bidirectional Sanger sequencing. We edited and assembled chromatograms and aligned sequences in Geneious 11.0.3.

(Biomatters Ltd).
First, we amplified the mitochondrial cytochrome oxidase I (COI) gene for all 864 specimens, for purposes of delimitation. We built a dataset of COI sequences of all known Niphargus species, regardless their origin , combined with all previously and de novo sequenced Western Balkans specimens from SubBioDB database, yielding COI sequences for 1492 individuals (1082 COI haplotypes) (Appendix S1). We used the alignment to delimit MOTUs, using distance-and tree-based unilocus delimitation methods: Automatic Barcode Gap Discovery (ABGD) (Puillandre et al., 2012), Assemble Species by Automatic Partitioning (ASAP) (Puillandre et al., 2021) and Poisson Tree Processes (PTP) (Zhang et al., 2013), respectively. In previous studies, we detected poor performance of multi-rate PTP (Delić et al., 2020) and discarded this method for Niphargus dataset.
Automatic Barcode Gap Discovery (Puillandre et al., 2012) assigns the sequences to the specific MOTUs without a priori species hypotheses, based on the assumption that intraspecific genetic distances are smaller than interspecific ones. We ran two ABGD analyses: with Kimura two-parameter substitution model, prior for maximum value of intraspecific divergence between 0.001 and 0.1, 20 recursive steps and two gap widths of 1.0 and 1.5. We used ABGD web server (https://bioin fo.mnhn.fr/abi/publi c/abgd/abgdw eb.html). For ABGD, we trimmed the alignment to 467 base pairs and excluded 16 sequences that were too short.
Similarly to ABGD, Assemble Species by Automatic Partitioning (ASAP) (Puillandre et al., 2021) is based on pairwise genetic distances, but provides a score for each defined partition and overcomes the challenge of a priori defining maximal genetic intraspecific divergence P. We ran ASAP with Kimura two-parameter substitution model, on ASAP web server (https://bioin fo.mnhn.fr/abi/publi c/ asap/asapw eb.html).
An alternative approach is a phylogeny-based method PTP (Zhang et al., 2013) that delimits species at nodes where presumed intraspecific nucleotide substitution rates switch to interspecific substitution rates. The two nucleotide substitution rates are modelled using two different Poisson processes. For this analysis, we removed duplicate sequences, using subset of 1082 unique haplotypes. First, we inferred phylogenetic relationships with maximum likelihood approach in IQ-TREE 1.6.7 (Minh et al., 2020;Trifinopoulos et al., 2016), using automatic best fit substitution model search with ModelFinder (Kalyaanamoorthy et al., 2017). We used the resulting tree to run the PTP analysis within the Bayesian and maximum likelihood framework, using the species delimitation server http://speci es.h-its.org/, and running 500,000 generations, sampling every 100 generations, and discarding the first 20% of the samples as a burn-in.
The delimitation results of the 1212 individuals from study area suggested between 212 and 474 MOTUs in the study area (see Results and Table 1). In order to reconcile the alternative results of the species delimitations for the downstream phylogenetic and spatial analyses, we proceeded as follows. We assigned five putative MOTUs statuses to each specimen. These results were crosschecked and compared with previous taxonomic studies. These studies employed multilocus species delimitations and/or additionally analysed morphological variation Delić et al., , 2022Fišer et al., 2015;Šet & Borko, 2020;Zakšek et al., 2019), and we deem taxonomic conclusions of these studies more reliable than unilocus delimitations alone. We first determined unproblematic MOTUs that were recognized by all five delimitation methods and, if applicable, agreed with previous works (478 out of 1212 specimens). Next, we revised the conflicting results. Whenever available, we chose delimitation results that corroborated with previous taxonomy-focused studies (433 out of 1212 specimens). If no such previous study existed, we assigned the specimens to a MOTU that was prevalent in the five species Final determination based on crosschecking with multilocus/integrative approaches 227 Note: ASAP appeared to perform the best (87% specimens delimited correctly, based on our criteria (see Methods section)). a The two ABGD analyses differ in pre-set gap width. We tested two settings that proved useful in different clades in previous analyses.
TA B L E 1 Delimitations based on COI marker, using ABGD, ASAP and PTP, reported for subset of Dinaric COI sequences F I G U R E 1 Workflow of the applied delimitation procedure. We combined existing and newly obtained COI sequences into an alignment. We applied five delimitation methods: ASAP, ABGD with two settings, mlPTP and bPTP, and assigned respective MOTU statuses to each individual. The results were compared with previous integrative taxonomic studies (Quality Check, QC: multilocus species delimitations and/or studies combining multiple molecular markers with morphology and spatial information). We first determined MOTUs that were recognized by all five delimitation methods unanimously (478 individuals). Next, we revised the individuals that were assigned to different MOTUs, and one assignment prevailed. 231 of them were confirmed with QC. For 256 QC was not possible, and majority rule prevailed. In the remaining conflicting cases without majority rule to follow, we either followed QC whenever possible (217)  delimitation methods (271 out of 1212 specimens) or to more conservative MOTU in cases where two MOTUs were equally likely (30 out of 1212 specimens). MOTUs that could not have been reconciled with previous multilocus delimitations were conservatively excluded from the dataset and are not accounted for in the summary report (35 specimens). Figure 1 summarizes the workflow of delimitation procedure. This procedure yielded 227 MOTUs.
This dataset was complemented with species for which we had DNA fragments other than COI (mostly 28S), but were genetically clearly distinct (i.e. independent lineages on the phylogeny). The final dataset was comprised of 245 MOTUs from the Dinaric Region, present in altogether 598 localities (see Appendix S2 in Supporting Information).
For purpose of phylogenetic inference, we selected specimens for additional sequencing in a way that each MOTU was covered with at least one specimen. For this subset, specimens representing MOTUs were amplified for seven additional markers: two segments of 28S rRNA gene, the histone H3 gene (H3), a part of the ITS rRNA gene, as well as partial sequences of the glutamyl-prolyl-tRNA synthetase gene (EPRS), heat shock protein 70 (HSP70) and arginine kinase (ArgKin) genes. Sequencing, editing and assemblage procedures were the same as described above. Oligonucleotide primers and amplification protocols were the same as in .

| Phylogeny and molecular clock
To assess phylogenetic diversity patterns, we built a multilocus timecalibrated phylogeny of all Niphargus MOTUs (see Appendix S3 in Supporting Information). For each marker, we aligned the sequences in Geneious, using MAFFT v7.388 plugin (Katoh & Standley, 2013), with E-INS-i algorithm with scoring matrix 1PAM/k=2 and the highest gap penalty. Alignments were concatenated in Geneious. To remove gap-rich regions from the alignments of non-coding markers, we used Gblocks (Talavera & Castresana, 2007) with less stringent selection set. To determine optimal substitution model for each partition, we used Partition Finder 2 (Guindon et al., 2010;Lanfear et al., 2012) under the corrected Akaike information criterion (AICc).

We built a time-calibrated chronogram of all known Niphargus
MOTUs with BEAST 2. We used the same four calibration points that proved to be congruent with each other (details in Borko et al., 2021).
Briefly, we assigned 1) a "modern-looking" Niphargus fossils from Eocene Baltic amber to the node where fossil's characters evolved for the first time; 2) the age of final submergence of the land bridges between Eurasia and North America to a root node given that Niphargus does not live in Nearctics; 3) the age of the last land-bridge between Crete and Greek mainland to the node in which monophyletic clade on Crete split from Greek taxa; 4) the opening of the connection between Paratethys and the Mediterranean basin to the node of the monophyletic Middle East clade, which derived from Eastern European taxa.
For each gene partition, we used a set of priors as followed: linked birth-death tree model, unlinked site models as in previous analyses, with fixed mean substitution rate and relaxed clock log-normal distribution with estimated clock rate. We used default settings of distributions of all estimates. We run the analyses for 200 million generations, sampled every 20,000 generation. The first 25% of trees were discarded as burn-in.

| Diversity patterns
To assess the groundwater biodiversity patterns within the Western Balkans, we explored taxonomic and phylogenetic diversity of Niphargus in space. We used a grid of 415 quadratic cells with 20 × 20 km resolution that has been evaluated as most suitable according to procedure in Zagmajster et al. (2008) (results not shown) and has been used in previous studies of subterranean diversity in the region (Bregović et al., 2019;Bregović & Zagmajster, 2016;Zagmajster et al., 2008), with the Lambert Conformal Conical Projection (central meridian 18°, parallels 42° and 46°).
Taxonomic diversity was calculated as sum of different MOTUs occurring within each grid cell, i.e. species richness (SR). We plotted SR on a raster of study area.
Phylogenetic diversity was computed using Faith's phylogenetic diversity (FPD), which is the sum of the branch lengths of a phylogenetic tree connecting all the species of a given assemblage. For calculation, the distance from the deepest connecting node to the root was excluded, resulting in single-species cells with zero FPD.
The alternative approach, which takes into account also branches of shared history, is more useful for phylograms which reflect individual species' accumulation of genetic changes, rather than chronograms where all species have same "time" distance to the root. We plotted FPD on a raster of study area.
Given that FPD positively correlates with number of species, we additionally analysed the pattern of standardized effect size of phylogenetic diversity (sesPD), a metrics of phylogenetic diversity that takes into account number of species per cell. sesPD compares empirical FPD against expected FPD that is estimated from a null model. Null model was generated in 9999 randomizations, in which species were randomly rearranged among the cells, while maintaining species' occurrence frequency and species richness, i.e. number of species per cell (independent swap algorithm (Gotelli, 2000)).
Phylogenetic diversity of a cell expected under the null model is a mean FPD calculated from 9999 randomly obtained FPDs. This procedure allows also an estimation of p-value, saying in which cells the observed FPD significantly deviate from expected FPD in the assemblage under null model of species distribution. Low sesPD is expected when two or more closely related species occupy the same cell, and might point to local speciation. High sesPD detects co-occurrence of distantly related species and may be pointer of dispersal causing cross-sections of clade's ranges. We calculated standardized effect size of phylogenetic diversity (sesPD = (observed PD -mean null PD)/standard deviation of null PD) and p-value as a quantile of observed PD vs. null communities (pd.obs.p = mpd.obs.rank/ runs +1). We plotted sesPD on a raster of study area and marked cells with significant p-value. Additionally, we applied another method to identify global variation in the PD controlling for the SR: For each cell with more than one species, we regressed PD against SR (see Appendix S4 for summary of linear regression) and calculated the residuals (the difference between observed and predicted PD, further on referred as residual PD). Cells harbouring more PD than expected for the observed richness are considered as regions of disproportionately phylogenetically diverse species compositions and vice versa (Gumbs et al., 2020). We plotted the residual PD on the 20 × 20 km grid covering the study area. Last, we compared the relationship between SR and sesPD. We plotted sesPD against SR and marked species-rich cells (seven or more species) based on their standardized phylogenetic diversity, to distinguish cells with lower and higher phylogenetic diversity than the mean value of null models, respectively. We plotted those cells on a raster of the study area.

| Taxonomic structure of the dataset and phylogenetic analyses
The results of five delimitation methods suggested that our Western Balkans dataset comprised between 212 and 474 molecular operational taxonomic units (hereafter MOTUs) ( Table 1). The difference in MOTU number appeared to be mostly due to different degrees of splitting, and in lesser extent due to incongruent lumping of different individuals into a single MOTU. After crosschecking with published multilocus species delimitations, we designed the dataset containing 227 MOTUs which can be considered as distinct species. After addition of 18 MOTUs without available COI marker, the final count was 245 MOTUs from Dinaric Region. ASAP outperformed other unilocus delimitation methods for large Niphargus dataset with 87% specimens classified to MOTUs that we consider as distinct species (Table 1).
To evaluate the molecular coverage of nominal species, we compared checklists of nominal species from that area (World Amphipoda Database (Horton et al., 2021)) with our data. Of the 123 nominal species listed from the area, we sequenced 97 (79%) species. Conversely, molecular delimitations indicate that the true species richness of the region may be 2 to 2.6 times higher.
A time-calibrated multilocus phylogeny of the whole genus included 512 MOTUs (Figure 2, Appendix S2). Of those, 377 were already recognized in previous studies  and 135 that were identified for the first time. Phylogeny comprised several geographically distinct, well-supported and MOTU-rich clades.
Most of the species in the study region belonged to nine of these clades that started to diversify approximately 15 million years ago ( Figure 2).

| Spatial analyses
To assess the biodiversity patterns in the study region, we mapped taxonomic and phylogenetic diversity of Niphargus using a grid of F I G U R E 2 Chronogram of Niphargus, obtained with BEAST2. Grey dots mark nodes with posterior support >0.9. Red dots on tips are MOTUs present in West Balkans. On the map in the middle localities of MOTUs used for phylogenetic inference are marked with either red (within West Balkans) or black dots (outside study area). Clades with more than five MOTUs from Western Balkans are highlighted, and maps at the outer side of the chronogram show their distributions inside the study area (study area is marked with dark grey)

| DISCUSS ION
The contemporary biodiversity crisis, with the number of endangered species far exceeding available conservation resources, requires knowledge of exceptionally biodiverse areas to maximize the impact of conservation acts (Myers et al., 2000). Here, using an exceptionally biodiverse system of subterranean amphipods, we assessed how historical processes may explain the discrepancy among biodiversity patterns of different metrics and how these processes could be important for conservation-oriented activities.
We identified different biodiversity patterns of Niphargus when comparing species richness and phylogenetic diversity within the region. The pattern of species richness, inferred from MOTUs, resembles the pattern reconstructed from the distributional data of morpho-species (Bregović et al., 2019). Although the datasets differ in overall number of data records and taxonomic units used for assessment of the biodiversity patterns, the analyses of both datasets agree that species richness peaks in the north-west, but few species-rich cells can also be found in the south-east. Noteworthy, similar patterns of higher and lower diversity in the north and south were detected also on the level of individual Niphargus clades (Delić, Švara, et al., 2017;Zakšek et al., 2019). Comparable species diversity patterns derived from different datasets of different hierarchic levels indicate that our results are solid, even though we acknowledge that both parts of the region might unveil additional, not yet detected species, or extend the geographic coverage to additional (yet empty) cells. It is noteworthy that past studies showed that hotspots in subterranean biodiversity remain hotspots even with additional sampling and that new species can be expected in both species-poor and species-rich areas (Culver et al., 2004;Zagmajster et al., 2008Zagmajster et al., , 2010. The pattern of phylogenetic diversity, however, is different: Phylogenetic diversity is low over most of the region except in the north-west, on the junction of the Dinaric Karst, South-Eastern Alps and Pannonian lowlands, at the crossroad of karstic areas and alluvial plains (Figure 3). The increase in phylogenetic diversity in the northwest can be attributed to the presence of species from two clades,

Pontian and Pannonic clades that are distributed between Northern
Italy and Pannonian lowlands, and cross Dinaric region in the north ( Figure 2, see also Borko et al., 2021). Europe-wide sampling  suggests that it is highly unlikely that members of either of these two clades would be found in the south-eastern part of the region. Likewise, we did not detect some clades from the south-east in the north-west; hence, both subregions contain some exclusive clades. We conclude that phylogenetic diversity results can be considered robust.
The patterns of species and phylogenetic richness seem to be decoupled. Spatial distribution of both, sesPD and residual PD, implies that species-rich and species-poor cells may have high PD (Figures   3c,d and 4). The contrasting differences among patterns of species richness and phylogenetic diversity of Niphargus in the Western Balkans allow inferences of the processes behind the pattern. Here we limit our discussion to cells with higher species richness in the north-western and south-eastern parts of the region. Species-rich areas can emerge either through dispersal from neighbouring regions or from on-site speciation (Rosindell & Phillimore, 2011). The pronounced roles of speciation and dispersal can be detected on patterns of phylogenetic diversity. Local species proliferation only little increases phylogenetic diversity and results in low phylogenetic diversity relative to species richness. Conversely, immigration of unrelated phylogenetic lineages increases phylogenetic diversity relative to species richness (Davies et al., 2007;Fritz & Rahbek, 2012;Li & Yue, 2020). Therefore, we hypothesize that speciation can be considered as a universal generator of species richness patterns along the entire Western Balkans and that dispersal acted mostly in the north-west (Figure 3). This view is consistent with the observation that distributional ranges of Niphargus species rarely exceed 200 km in length (Bregović et al., 2019;Trontelj et al., 2012). We recognize that the differences in species and phylogenetic richness could be also an outcome of extinctions. While we cannot directly evaluate the latter's contribution with the data available, we assume that this mechanism was negligible and random in space, given that subterranean environment acts as a refugium from major catastrophic disturbances on the surface . The most recent driver of extinctions in Europe were Pleistocene glaciations, but at that time the entire Balkan Peninsula acted as a southern refugium (Hewitt, 2000). Even more, Alpine glaciers in the north-west of the studied area prompted Niphargus speciation rather than extinction (Delić et al., 2022). The Balkan Peninsula has a dynamic geological history (Handy et al., 2015;Popov et al., 2004), which directly shaped the evolutionary trajectory of Niphargus in the area. An uplift of carbonate mountain ranges-including the South-Eastern Alps and the Dinaric Karst (Park, 2014;Popov et al., 2004)-opened up a variety of new, unpopulated habitats, acting as an ecological opportunity for rapid ecological diversification of several Niphargus clades . Highly fragmented subterranean environment also generated dispersal barriers that might have promoted allopatric speciation, resulting in cryptic species Zakšek et al., 2019). These adaptive and non-adaptive local radiations resulted in cells with low standardized phylogenetic diversity, especially in the south-east. Nevertheless, some Niphargus species from non-Dinaric clades have large ranges and are good dispersers (Copilaş-Ciocianu et al., 2018). Those species migrated in the north-western Dinarides . The north-western part of the studied region was occasionally submerged (Kováč et al., 2018). Aquatic connections be- species-rich areas may be an assembly of sink populations rather than regions with high functional diversity that stabilizes ecosystem functioning (Funk & Burns, 2019). If so, proper conservation strategy should not protect only the species-rich focal region, but also contributing areas.
A completely different issue is found for cells with low phylogenetic diversity, for which we hypothesize that gained species richness mainly by speciation, e.g. through local adaptive radiations . Such cells with high species richness and low phylogenetic diversity can be found in the south-east (Figure 4).
Anthropogenic disturbance in these areas may negatively affect whole lineages, ranging from loss of diversity by reverse speciation (Seehausen, 2006), to extinction of entire clades and irreversible loss of their genetic diversity. Such locally distributed and closely related species complexes can be more threatened than unrelated species complexes simply because they are confined to a threatened environment (Funk & Burns, 2019;Tonini et al., 2016) or also because of inherited sensitivity (phylogenetic niche conservatism) (Wiens et al., 2010). Potential threats could damage the habitat template that had prompted diversification in the past, and thus make recovery of the local ecosystems less likely. plants present a direct threat to the enigmatic biodiversity of this region and Europe as a whole (Schwarz, 2015).
Identification and comprehensive analysis of contrasting biodiversity patterns sheds a new light on the origin of species-rich areas of the Western Balkans. We here identify several issues that should be addressed when planning conservation measures: the hidden biodiversity (undescribed, sometimes morphologically cryptic species); the discrepancy between species richness and phylogenetic diversity patterns; and possible implications of different origins of hotspots on conservation planning. There are some additional aspects that may be relevant for conservation biology.
Despite massive efforts in the past decades, we can expect additional species records (Appendix S5, Figure S5.1) and undescribed species (Table 1). Second, although Niphargus is an important faunistic element of the studied region, it is merely one among many taxa, many of which are still understudied. It would be beneficial to include other taxa, aquatic and terrestrial, which may reveal different relationships between species richness patterns and phylogenetic diversity. However, identification of robust biodiversity patterns is only the first step, preceding a challenging task how such science backed knowledge should be integrated into overarching conservation strategy in the region spanning across four countries. Finally, this study implies that similar approaches should be applied also in other regions, at least where species richness is high (Eme et al., 2018;Zagmajster et al., 2014). However, the generally scarce sampling of groundwater, with only few studies having similarly dense spatial resolution (Alther et al., 2021;Fišer et al., 2018;Flot et al., 2014;McInerney et al., 2014), makes such inclusions unlikely to happen in the near future. Thus, while we are still far from an optimal strategy on how to protect groundwater habitats, the numerous threats to groundwater (Mammola et al., 2019), such as hydropower dams or extensive water extraction, put an extra urge to strengthen and systematically organize future research on this topic, to make it implementable and effective. and Biodiversity. We thank four anonymous reviewers for their comments.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sequence data have been deposited in GenBank, under accession numbers OK149784-OK150101 and OK156503-OK157282. Author Contribution: Š. B., M. Z., F. A. and C. F. conceptualized the study. Š. B. and M. Z. did the field work. Š. B. performed laboratory work and analysed the data. Š. B. prepared the first draft.
Š. B., M. Z., F. A. and C. F. wrote and approved the paper.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found in the online version of the article at the publisher's website.