Can Environment Predict Cryptic Diversity? The Case of Niphargus Inhabiting Western Carpathian Groundwater

In the last decade, several studies have shown that subterranean aquatic habitats harbor cryptic species with restricted geographic ranges, frequently occurring as isolated populations. Previous studies on aquatic subterranean species have implied that habitat heterogeneity can promote speciation and that speciation events can be predicted from species’ distributions. We tested the prediction that species distributed across different drainage systems and karst sectors comprise sets of distinct species. Amphipods from the genus Niphargus from 11 caves distributed along the Western Carpathians (Romania) were investigated using three independent molecular markers (COI, H3 and 28S). The results showed that: 1) the studied populations belong to eight different species that derive from two phylogenetically unrelated Niphargus clades; 2) narrow endemic species in fact comprise complexes of morphologically similar species that are indistinguishable without using a molecular approach. The concept of monophyly, concordance between mitochondrial and nuclear DNA, and the value of patristic distances were used as species delimitation criteria. The concept of cryptic species is discussed within the framework of the present work and the contribution of these species to regional biodiversity is also addressed.


Introduction
In the last decade, the rise of molecular studies has greatly improved the detection of cryptic species (morphologically indistinguishable species) [1][2][3][4][5][6][7]. The importance of cryptic species lies in their contribution to overall biodiversity by increasing species richness at different scales and may also be due to the fact that they convey important information as fundamental units in biogeography, ecology, evolutionary studies [8], [9] and conservation biology [10][11][12]. Uncovering cryptic diversity is important for understanding species distribution ranges, assessing levels of endemism and for species ecology, as well as for the conservation status of such cryptic species [13], [14].
Some authors have suggested that the phenomenon of crypsis is rather ubiquitous across all animal phyla and regions [15]. On the other hand, a study firmly grounded in evolutionary theory has suggested that groups living in environments with strong directional selection might be subject to morphological crypsis more often [16]. The subterranean realm is a highly fragmented environment where strong directional selection operates. The fragmented nature of its habitats increases the possibility of speciation, whereas strong directional selection constrains the extent of morphological changes [5][6][7], [17][18][19][20][21][22][23][24][25][26][27][28]. Among aquatic subterranean taxa, cryptic species seem to be common [29].
A review of the published data suggests that breaks in gene flow can be inferred from (i) a geologically heterogeneous environment (N. rhenorhodanensis [17]), (ii) breaks among water catchments (Niphargus virei [21]; Proteus anguinus [30]; Troglocaris anophthalmus [18] or (iii) other types of environmental heterogeneity (N. ictus [27], N. rhenorhodanensis [31], N. virei [21]). Therefore, if the ubiquity of cryptic species remains elusive, can we at least predict which morphospecies more likely to consist of two or more cryptic species using environmental cues and species distributions?
We approach this issue using the case of the subterranean West Palearctic genus Niphargus. With more than 300 species, Niphargus is the most speciose freshwater amphipod genus in the world [32]. A high level of cryptic diversity within Niphargus has been uncovered by molecular studies [5], [21], [27], [31], [33]. However, there are no molecular studies for the Eastern part of Europe.
In this study we test a bold prediction, i.e. that cryptic diversity is to some extent predictable. In other words, morphospecies from ecologically heterogeneous environments and/or distributed across different drainage systems likely represent complexes of morphologically similar species. Our study took place in an important small conservation area that (i) has been thoroughly studied using traditional taxonomic approaches and that (ii) is heterogeneous and highly fragmented and therefore satisfies the requirements for environmental heterogeneity. The prediction of cryptic diversity was addressed by discussing trends in niphargid speciation in the Western Carpathians using a molecular approach on Niphargus populations from 11 caves distributed across the mountain range.

Sampling Area
One of the main characteristics of the Western Carpathians (the so-called Apuseni Mountains) is the relatively high percentage of karst landscape (11%) compared to other karst areas in Romania, which covers a surface of about 10750 km 2 , with an average elevation of 700 meters. The Apuseni Mountains were chosen as a study area due to the following characteristics: i) the diversity of Niphargus species based on morphology was well studied within this area, ii) the highly fragmented karst landscape [34] is prone to harbour a high level of cryptic diversity and iii) waters from the area have outflows into two different drainage systems (the Crişul Repede and Crişul Negru basins). We collected samples in 11 caves, focusing on percolation waters and pools. The sampled localities are shown in Figure 1.

Taxa Sampling
An overview of the published data (review in [35]) suggests that almost 20% of all Romanian niphargids are found in the Western Carpathians, representing almost half of all Niphargus species known from Romanian caves. A survey of the available literature revealed the presence of Niphargus species in 18 caves in the Western Carpathians. We sampled nine of those caves plus four  Table 1 Collected individuals were identified using identification keys and original species descriptions [36]. Altogether, 59 specimens were analyzed molecularly.

Molecular Protocols
Genomic DNA was extracted from one pereiopod (the rest of the animal was kept for morphometric studies) or from the whole specimen for small individuals using the GenElute Mammalian Genomic DNA miniprep kit (Sigma-Aldrich). The first 28S rDNA fragment was amplified using the forward primer from [37] and the reversed primer from [38]. For the second part of the 28S rDNA fragment, the pair of primers from [39] was used. The H3 histone was amplified using primers H3NF and H3NR from [40]. In addition, mitochondrial cytochrome oxidase I (COI) was amplified using LCO1490 and HCO2198 primers [41], which is a widely used barcoding marker for genetic diversity within and between closely related populations.
PCR was performed using the following cycling settings: 45 s at 94uC, 30 s at 48uC, 90 s at 72uC, for 35 cycles followed by a final extension at 72uC for 3 min (the first part of 28S rDNA); 45 s at 94uC, 60 s at 48uC, 120 s at 72uC, for 30 cycles followed by a final extension at 72uC for 3 min (the second part of 28S rDNA); 45 s at 94uC, 60 s at 46uC, 60 s at 72uC, for 35 cycles followed by final extension at 72uC for 3 min (H3); 60 s at 94uC, 60 s at 45uC, 150 s at 72uC for 40 cycles followed by a final extension at 72uC for 7 min (COI). PCR products were purified using Exonuclease I and Alkaline Phosphatase (Fermentas Inc., Germany). Each fragment was sequenced in both directions using PCR amplification primers from Macrogen Europe (Amsterdam, The Netherlands). Contings were assembled and edited using Chromas Pro Version 1.5 (Technelysium Pty Ltd).

Phylogenetic Analyses
We compiled four different datasets. The phylogenetic position of Romanian niphargids was inferred from the Niphargus dataset of the first 28S gene fragment and from the outgroup available in GenBank [5], [39], [42]. Details of the samples are shown in Table 1 and Table S1. To determine the speciation of specimens from the Western Carpathians we used the following datasets: COI, the second part of 28S, and H3, and combinations of these datasets (Dataset S1).
The sequences of different Niphargus species from the first part of the 28S rDNA varied considerably in length; the length of the 28S alignment was about 1100 bp. The large differences were mainly due to simple sequence repeat insertions in some species. To account for the long indels, the sequences of 28S were aligned using the E-INS-i option for sequences with multiple conserved domains and long gaps in MAFFT ver. 6 [43]. Low homology regions with long gaps were removed using Gblocks [44] under the least restrictive settings possible. Altogether, 1075 nucleotides were kept for phylogenetic analyses.
A general time-reversible model with a proportion of invariant sites and a gamma distribution of rate heterogeneity (GTR+I+C) assuming six discrete gamma categories was chosen as the most appropriate model according to AIC and BIC criteria, using ModelGenerator [45]. Bayesian analyses were performed using MrBayes v3.1.2. [46]. Two parallel searches with four chains each were run for two million generations sampled every 100 th generation. The burn-in value was graphically determined from the plot of the likelihood values of the trees. The trees visited by the chains before the likelihood values reached a plateau were discarded as burn-in. The final topologies were constructed according to the 50% majority rule. The maximum likelihood (ML) phylogeny was obtained using PHYML [47]. All parameters of the nucleotide substitution model and the gamma shape parameter were simultaneously estimated during the ML search. The robustness of the topology was tested with 1000 bootstrap support. Maximum Parsimony (MP) was performed in PAUP version 4.0b10 [48]. The best tree was searched using the heuristic algorithm, by performing tree bisection and reconnection with taxa added randomly with 10 replications. The robustness of the nodes was estimated by performing 100 bootstrap replicates for the large 28S dataset and 1000 bootstrap replicates for the Dataset S1. Phylogenetic analyses (i.e. ML) reconstruct evolutionary relationships of sequence data under the assumption that a tree represents their best relationship. For similar sequences, due to low variability at the intraspecific level, these relationships are often more clearly and accurately represented by networks [49]. Therefore a median-joining network [50] was conducted within the ''Laticaudatus'' clade on the COI dataset (33 specimens from five populations). The analysis was performed using the program Network ver. 4.6, (available at www.fluxus-engineering.com). We assigned equal weights to all positions and e was set to zero.

Genetic Divergence
Molecular divergences were calculated using patristic distances. Patristic distances were calculated using the PATRISTIC v1.0 program from an ML tree as described in [51].

Species Delimitation Criteria
An accurate delimitation of species is essential as species are the basic units for biodiversity studies and as such they are the basic units for conservation strategies. We have defined species as independently evolving lineages, i.e. a general species concept introduced by de Queiroz [52], [53]. In time, these lineages may evolve a range of characteristics like morphological and genetic distinctness that result in exclusive monophyly, sometimes in ecological distinctness, and finally in the evolution of a reproductive barrier. The species delimitation criteria that we have used include the concept of monophyly, concordance between mitochondrial and nuclear DNA, and the threshold value of patristic distances (for details on the criteria see publications [54], [55]). We applied these criteria in order to assess the evidence for or against the species status of an individual population/group of populations.

Morphological Identification and Selection of Characters for the Illustration of Morphological Similarity
Firstly, we identified specimens according to the available descriptions and diagnoses [36]. Using this information, we diagnosed four distinct species: Niphargus andropus, N. laticaudatus, N. transsylvanicus and N. bihorensis, although we were not able to collect N. stygocharis and N. stygius, which are known to be present in Vadu Crişului Cave and Ferice Cave [56], [57], respectively. N. andropus is a small species that was represented in only a few samples by one or two individuals. Due to their small size, whole individuals were used for DNA extraction, therefore later morphological examination was not possible. Similarly, N. transsylvanicus was represented by a single individual. In the remaining two species (N. laticaudatus and N. bihorensis) where, after molecular analyses, it turned out that they comprise a complex of species, we searched for additional morphological differences between them. We checked for 10 continuous morphological characters commonly used in Niphargus taxonomy, including those that turned out to be useful for species delimitation in morphologically similar species [58][59][60][61].
Continuous characters included measures of the head and first pereonite (a surrogate for the body length), length of antenna I, six measures on gnathopod II (depth and width of coxal plate, length of carpus, length of propodus, palm length of propodus, distance between palmar spine and carpo-propodal article) and two measures on pereopod VII (length of appendage and width of article 2). Details about characters and landmarks are presented and discussed in [62].
For N. bihorensis species complex, we analyzed 10 adult individuals from each population. In N. laticaudatus we used adult individuals from Corbasca Cave (six individuals), Ferice Cave (six individuals), Grueţ Cave (five individuals), Osoi Cave (six individuals) and Ungurului Cave (six individuals). The characters we used are sexually non-dimorphic, therefore we included both sexes.
Niphargus specimens were partially dissected in glycerin and mounted on slides. Appendages were photographed with an Olympus camera ColorView III mounted on an Olympus DP Soft stereomicroscope and measured using the ANALYSIS (Olympus Soft Imaging Solutions) program. The appendages and the rest of body were stored in the Zoological Collection of the Department of Biology, Biotechnical Faculty, University of Ljubljana.

Statistical Analyses used for Morphometric Data
Continuous characters were analyzed together, using principal component analysis (PCA) on a covariance matrix. Some specimens were partially damaged. Missing values were replaced Table 2. Species delimitation criteria met by niphargid species (according to [54], [55] using the expected value for the trait as estimated from the regression line (trait-body size) calculated from conspecifics. The first two principal components were plotted for visual inspection if pairs of cryptic species showed any grouping in morpho-space. The analysis was performed using PASW ver.18 software.

Phylogenetic Relationships
Maximum likelihood, maximum parsimony and Bayesian inference were conducted on the dataset of the first 28S gene fragment of 87 Niphargus species and resulted in trees with similar  Table 1 Figure 2). The phylogenetic analyses of Niphargus species from the Western Carpathians within the genus Niphargus showed that they belong to two completely different and phylogenetically independent clades. We designated them as ''Andropus-Bihorensis'' and ''Laticaudatus'' clades ( Figure 2). The ''Andropus-Bihorensis'' clade included three morphospecies from seven caves: Ciur Izbuc, cu Apȃ din Valea Leşului, Drȃcoaia, Mȃgura, Meziad, Osoi and Vadu Crişului. This clade comprises species distributed from from Western Europe (France, Northern Italy) to Central and Eastern Europe. The area between Northern Italy and Romania has been invaded by species from two subclades, both of them reaching the Western Carpathians. For clarity, we have designated them as ''Bihorensis'' and ''Andropus'' subclades. The two Carpathian subclades are more closely related to other non-Carpathians species than to each other. The ''Bihorensis'' subclade contains N. bihorensis whereas ''Andropus'' contains N. andropus and N. transsylvanicus (Figures 2, 3, 4, Table 2). Both subclades are genetically more diverse than morphologically and in fact comprise two and four species, respectively (see below).
The other completely separate clade of Niphargus from the Western Carpathians that we called ''Laticaudatus'' was taxonomically much less diverse than ''Andropus-Bihorensis'' clade and consists of a single morphospecies distributed across five caves. The relationship of the ''Laticaudatus'' clade to the rest of the Niphargus species is not certain (Figure 2). The ''Laticaudatus'' clade is also genetically more diverse than morphologically ( Figure 5, Table 2); northern populations (Ungurului and Osoi caves) have turned out to be genetically distinct from the southern ones (Ferice, Grueţ and Corbasca caves) (see below).

Genetic Divergence
A total of 668 nucleotides of the COI fragment were obtained from 55 Niphargus specimens from nine different localities in the Western Carpathians. The sequencing of COI for two populations unfortunately failed (see Table 1).
Patristic distances between the two pairs of species are as follows: the distance N. bihorensis -Niphargus sp. 4 is 0.19; the distance Niphargus sp. 2 -Niphargus sp. 3 is 0.21. The lowest patristic distance on COI is 0.04 between northern and southern ''Laticaudatus'' ( Table 2). All species for which both 28S fragments and H3 were analyzed have unique 28S (both fragments) and H3 sequences, satisfying the criteria of exclusivity and congruence between the independent markers ( Table 2). The taxonomic conclusions are summarized in Table 3.
In order to confirm or refute the existence of an additional cryptic species within the ''Laticaudatus'' clade, due to a relatively small variability on 28S and a low patristic distance between the northern and southern group, a detailed haplotype analysis was performed for the ''Laticaudatus'' clade. The haplotype network for ''Laticaudatus'' showed the split between the northern and southern populations as already indicated by nuclear data (Figure 2). The network analysis resulted in nine median vectors and three cycles ( Figure 5). It uncovered two groups of haplotypes separated by 22 mutational steps. As such, this indicates that the northern and southern populations are genetically separated and there is no indication of gene flow between them, as we had already predicted from the phylogenetic tree based on 28S ( Figure 2). Finally, we checked the sequences of an additional 28S fragment, which consistently supports the split between northern and southern populations (Dataset S1).

Morphological Analyses
All measurements and observations are found in Table S2. The first principal component (PC) chiefly explains body size variation (Table 4), while the second PC accounts for variation in appendage length and body shape. Both PCs together explain 97.80% and 99.40% of the total variation in the N. bihorensis and N. laticaudatus species complex, respectively. Examination of the plots shows that variation in both species pairs largely overlaps ( Figure 6).

Speciation in Niphargus from the Western Carpathiansthe Role of Hydrography
The possibility of dispersal within highly fragmented limestone patches and aquifers along the Western Carpathians is limited, which could cause the isolation of populations, followed by genetic differentiation and speciation. The morphospecies N. bihorensis and N. laticaudatus appear to be widespread in both hydrographic basins and in both massifs within the Western Carpathians (Bihor and Pȃdurea Craiului). However, molecular tools suggest that in fact each hydrographic basin (Crişul Negru and Crişul Repede) is hosting cryptic sister species. Populations in each hydrographic basin are characterized by monophyly and concordance between nuclear and mitochondrial data. The patristic distance on COI between northern and southern ''Laticaudatus'' is below the proposed cut-off divergence level for species delineation in crustaceans [54], however it is within the cut-off divergence level proposed for marine amphipods [63]. Whatever the lower boundary cut-off divergence is, genetic and distributional data suggest that the two lineages within the ''Laticaudatus'' clade have been evolving independently for some time. We find the assumption that on-going speciation has been mediated by the hydrographic regime to be reasonable. The inference of speciation based on hydrographically separated areas seems to be straightforward; however, it is surprising that none of the N. bihorensis species can disperse through tiny voids in carbonate rock over such a short distance (20 km). One can infer that their genetic structure might be spatially associated with the Plio-Quaternary paleo-drainage network, which had two main paleo-flow directions in Apuseni: a North-East drainage corresponding to the more recent Crişul Repede hydrographic basin and a South-West drainage corresponding to the more recent Crişul Negru hydrographic basin [64]. The overall morphological similarity of the sister species suggests that the evolution of morphology is constrained by strong environmental selection.

Speciation in Niphargus from the Western Carpathiansthe Role of Fragmented Karst
The diversity of the ''Andropus'' subclade indicates that speciation may also occur within the same basin, which is consistent with previous studies on groundwater species [5]. Three different species in this subclade are located at most 30 km away from each other, all within the same Crişul Negru hydrographic basin, while the forth species was found in the Crişul Repede hydrographic basin. The taxonomic conclusions for the ''Andropus'' subclade remain tentative and we propose that it consists of four phylogenetic species: N. andropus, Niphargus sp. 2, Niphargus sp. 3 and N. transsylvanicus. The low number of analyzed individuals and the failure of DNA sequencing (on COI) in two out of three individuals made the classification somewhat difficult. However, four populations (Ciur Izbuc, cu Apȃ din Valea Leşului, Drȃcoaia and Mȃgura caves) were analyzed for 28S and histone H3 and all four are genetically distinct. The present distribution of the four populations can be explained by paleohydrographic processes that caused fragmentation by non-karstic deposits acting as natural barriers to species migration or by tectonic movements that changed the subterranean drainages [34], [65]. An alternative scenario for the within catchment speciation through isolation is the putative difference in ecology between Niphargus sp. 2 and Niphargus sp. 3. Niphargus sp. 2 exhibits small body size and was sampled only from dripping water, whereas the larger body size Niphargus sp. 3 was found in a large cave pool. This scenario is also supported by the small N. transsylvanicus (Crişul Repede hydrographic basin) sampled from a pool fed by percolating water, suggesting its habitat preferences for the fissure system within the limestone maze. This is in accordance with the observation that small-bodied species live in tiny crevices [66].
Distinct phylogeographic patterns of subterranean taxa at different taxonomic levels related to habitat fragmentation and heterogeneity are frequently mentioned in other studies of subterranean species, both aquatic and terrestrial [28]. In the Western Carpathians, cave-dwelling beetles are the only group of subterranean animals that have been molecularly analyzed to date [67]. That study based on the mitochondrial DNA of three genera has revealed that phylogeographic breaks within genera might be the result of karst fragmentation acting as long-term barriers to gene flow among species and subspecies. These findings are in accordance with the present data suggesting that using karst fragmentation to predict speciation patterns can be generalized for both terrestrial and aquatic species across the Western Carpathians.

Overlooked but Expected Diversity in the Western Carpathians -Implications for Biodiversity Research
This is the first study to reveal niphargid cryptic speciation in the Carpathians.
The obtained results are in agreement with our prediction that environment can predict speciation events, even when the morphology of distinct populations is not different. Both geological fragmentation as well as heterogeneous hydrogeological settings likely limit dispersal and promote speciation. Moreover, speciation should be expected within the entire genus since this is not a property of a certain clade with a strongly conserved ecological niche [68]. Partially predictable speciation has important consequences for biodiversity research. Using the available environmental layers and morphospecies distributions we can easily identify hypothetical morphospecies, i.e. species in which genetic diversity is most likely underestimated. In conclusion, linking morphospecies ecology with its distribution patterns could be used as a tool to reveal a more accurate picture of biodiversity across spatial scales. This would optimize taxonomic research with a relatively high degree of certainty, making taxonomy faster and more rewarding.

Supporting Information
Table S1 List of Niphargus taxa used in the phylogenetic analysis with their geographic origin and GenBank Accession numbers. (DOC) Dataset S1 Dataset S1. Phylogenetic trees conducted by three different phylogenetic methods: Bayesian inference (BI), maximum likelihood (ML) and maximum parsimony (MP) from different molecular markers and their combination of Niphargus from Romania. Bootstrap value (ML), posterior probabilities (BI) and bootstrap value (MP) are shown on each branch. Analyses were performed as described in the Materials and Methods section of the manuscript. (DOC)