Range sizes of groundwater amphipods (Crustacea) are not smaller than range sizes of surface amphipods: a case study from Iran

The connectivity of groundwater aquifers is lower compared to surface waters. Consequently, ground-water species are expected to have smaller distributional ranges than their surface relatives. Molecular taxonomy, however, unveiled that many species comprise complexes of morphologically cryptic species, with geographically restricted distributional ranges in subterranean as well as in surface waters. Hence, the range sizes of surface and groundwater species might be more similar in size than hitherto thought. We tested this hypothesis by comparing the range size of surface amphipods of the genus Gammarus and subterranean amphipods of the genus Niphargus in Iran. We re-analyzed the taxonomic structure of both genera using two unilocus species delimitation methods applied to a fragment of the COI mitochondrial marker, to identify molecular operational taxonomic units (MOTUs), and assessed the


Introduction
The groundwater is the largest reservoir of unfrozen freshwater in the World (Gibert & Deharveng, 2002).Despite that, groundwater metazoans count less than 10% of the global metazoan species richness (Gibert & Deharveng, 2002), suggesting that specific ecological conditions such as a permanent darkness, limiting food and reduced short-term climatic fluctuations (Culver & Pipan, 2009), act as a strong abiotic filter.Few taxa nevertheless successfully colonized and diversified in the groundwater (Trontelj, 2018), and nowadays represent some of the ecologically most specialized freshwater organisms and some of the oldest freshwater inhabitants (Humphreys, 2000;Culver & Pipan, 2009;Protas & Jeffery, 2012;Mammola, 2019).
Subterranean aquatic species have generally narrow distributional ranges, presumably due to the high degree of fragmentation of groundwater aquifers (Trontelj et al., 2009;Zagmajster et al., 2014;Delić et al., 2017;Eme et al., 2017).Species with large ranges are an exception (Eme et al., 2014;Copilaș-Ciocianu et al., 2017), and seemingly more common in northern latitudes (Zagmajster et al., 2014;Eme et al., 2017).Such narrow range-endemic species are generally considered as an important fraction of the global natural heritage.Their presence may steer nature-conservation decisions on a local and regional scale (Sket, 1999;Chapman et al., 2009).In comparison, surface habitats are better connected, the dispersal of surface species is less limited and these species have presumably larger distributional ranges (Cardoso, 2012;Altermatt et al., 2014;Alther et al., 2016).The differences in the within-surface and within-subterranean habitat connectivity and putative differences in range sizes imply a hypothesis that the subterranean fauna is on average more endangered and should be monitored more carefully than the surface fauna.
The hypothesis that subterranean species are geographically more restricted (endemic) and potentially more threatened than surface species, so far received little attention (but see Cardoso, 2012;Mammola et al., 2019), and remains to be tested for aquatic fauna.Over the past decades it has become clear that a substantial fraction of species diversity can be identified only with the help of molecular methods (Bickford et al., 2007;Fišer et al., 2018).Morphologically identical or highly similar species, so-called cryptic species, are common in all animal phyla and can be expected in all regions of the world (Pfenninger & Schwenk, 2007;Pérez-Ponce de León & Poulin, 2016).Many studies unveiled that also geographically widespread surface species in fact comprise complexes of narrow endemic species (Copilaş-Ciocianu & Petrusek, 2015;Mamos et al., 2016), thereby suggesting that the degree of endemism in surface waters may be underestimated and that the difference in endemism between surface and groundwater animals is less distinct than previously thought.
We assembled genetic and distributional data on genera Niphargus and Gammarus from Iran.Several samples of Niphargus from Azerbaijan province in Iran have been collected for the first time.To test our hypothesis, we applied genetic species delimitation methods to standardize taxonomic units, and defined range sizes as the maximum linear extent (MLE) between the two distant-most points.

The data collection
The samples were derived from various collections.Data for the genus Gammarus (supplementary table S1) were extracted from the database published in Katouzian et al. (2016).
The total dataset comprised 67 localities in the Irano-Anatolian and Caucasus area, and 180 COI sequences which yielded 104 haplotypes.

Phylogenetic analyses of the new samples
We extracted the total genomic DNA from pleon using Tissue Kits (GenNet Bio™) following the manufacturer's instructions (Seoul, South Korea).Mitochondrial COI was amplified using the modified primer pair LCO1490-JJ and HCO2198-JJ (Astrin & Stüben, 2008).Amplification and sequencing of the first fragment of 28S ribosomal DNA (rDNA) were performed using the forward primer developed by Verovnik et al. (2005) and the reverse primer used by Zakšek et al. (2007).Each 25 μl reaction consisted of optimized amounts of PCR water, 12.5 μl of Master Mix kit (Sinaclon, Iran), 0.2 μl of each primer (10 μM), and 50-100 ng of genomic DNA template.For COI gene amplification, an initial denaturation step at 94°C for 3 minutes was followed by 36 cycles of 40 seconds at 94°C, 40 seconds at 52.5°C and 2 minutes at 65°C with a final extension step for 8 minutes at 65°C.Cycling parameters for the 28S rDNA gene were as follows: initial denaturation of 94°C for 7 minutes, 35 subsequent cycles of 94°C for 45 seconds, 55°C for 30 seconds, 72°C for 1 minute, and a final extension of 72°C for 7 minutes.Purification of PCR products and sequencing Downloaded from Brill.com09/17/2023 10:29:08AM via free access were commercially performed by Macrogen Inc. (Korea).Sequencing was performed with both primers mentioned above.
Phylogenetic reconstruction was performed using Bayesian inference in Mr Bayes, version 3.1.2(Ronquist & Huelsenbeck, 2003).Analyses were run for five million generations, with four chains, and the trees were sampled every 1000 generations under GTR + G and GTR + I + G models (jModelTest, version 0.1.1,Posada 2008) for 28S and COI genes, respectively.The first 1250 sampled trees were discarded as burn-in, and the subsequent tree likelihoods were checked for convergence in Tracer 1.5.0 (Rambaut & Drummond, 2009).A fifty percent majority rule consensus tree was computed using the remaining trees and visualized by FigTree v1.4.0 software.The data on analyzed species are available in supplementary tables S1 and S2.

Species delimitation methods
We applied two molecular species delimitation approaches, namely Automatic Barcode Gap Discovery (ABGD), Bayesian Poisson Tree Process (bPTP) methods (Puillandre et al., 2012;Zhang et al., 2013).The same methods and settings were used as in previous studies of genus Gammarus populations (Katouzian et al., 2016), which we re-used in the analyses of genus Niphargus populations as follows.
The ABGD delimitation was run using the following settings: the parameter range of P min = 0.001, P max = 0.10, gap width = 0.8-1.1, and a Kimura-2-parameter (K2P) model (Kimura, 1980) corrected genetic distance matrix generated in MEGA ver. 5 (Tamura et al., 2011).In bPTP species delimitation were performed using the webserver (http://www .species.h-its.org)and default settings.As input, we used a COI tree produced with BEAST v1.7.4.The input file for BEAST was produced in BEAUTi v.1.7.The settings were as follows: 50 million generations, sampling each 5000th tree, standard coalescent model, GTR + I + G substitution model with four gamma categories and a strict clock.GTR + I + G model was proposed by model selection in jModel test.Appropriateness of parameters (effective sample size >200) was tested with Tracer v1.6.Results were visualized with Tree Annotator v1.7.4 with a 10% burn-in rate, posterior probability of 0.9 and under the maximum clade credibility option for the consensus tree.
For five Niphargus species we had only 28S sequences, and previous studies proved to be genetically and morphologically distinct from other Niphargus species (Hekmatara et al., 2013;Esmaeili-Rineh et al., 2016, 2017a, 2018).These species, however, are known from single sites and could be unambiguously included into the dataset.

Distributions
The maximum linear extents (MLE) for molecularly delimited operational units were obtained in ArcGIS (version 10.2).Species, known from a single locality were arbitrarily assigned a distributional range of 1 km.The distribution of MLE was left skewed and not normal, hence we tested whether ranges of surface and subterranean species differed in sizes using the nonparametric Manny-Whitney U test using IBM SPSS 20.
The genus Niphargus counted 16 described species and several putative new species which are pending further taxonomic work.New samples from west Azerbaijan province comprised two strongly supported subclades, both nested in the main Iranian clade (fig.2).The ABGD and bPTP species delimitations suggested that the samples collected from west Azerbaijan province comprised 8 and 9 species, respectively (fig.3; supplementary table S2).All known samples of genus Niphargus counted 24 and 26 MOTUs, according to ABGD and bPTP delimitations, respectively.The MLEs of genus Niphargus MOTUs ranges measured between 1 and 644 km and between 1 and 844 km, when we used bPTP and ABGD delimitations, respectively.

Discussion
The results of our study did not support the hypothesis that subterranean species have smaller distributional ranges than surface species and that they should be treated as more endangered per se.While the connectivity of surface amphipod populations within a single catchment is not problematic, the connectivity between catchments seems to be as challenging as the connectivity between groundwater aquifers.Amphipods carry their eggs in the marsupium and do not have desiccation resistant propagules (e.g., Kaestner, 1970), which might explain their relatively limited range sizes among surface species.By contrast, the groundwater aquifers may temporarily communicate through occasional subterranean connections, occurring at high water levels that at least temporary enhance the subterranean connectivity (Palandačić et al., 2012;Konec et al., 2016).Despite that small ranges are the rule in groundwater (Trontelj et al., 2009), in continuous karst areas, some species apparently can establish and maintain large distributional ranges using subterranean connections that do not reflect hydrological conditions on the surface (Lefebure et al., 2006;Delić et al., 2017;Zakšek et al., 2018).
Our results suggest that range sizes of surface and subterranean amphipod species should be revised in general.This conclusion is in accord with a recent review of cryptic diversity of the amphipod crustaceans (Fišer et al., 2018).The meta-analysis of studies suggested that cryptic amphipod species were found in all studied families, and from all environmental conditions, including surface and subterranean freshwater, semi-terrestrial, marine pelagial and benthos, deep sea and parasites.Virtually all these studies showed that distributional ranges of molecularly delimited species are much smaller than distributional ranges of morphologically defined species (Fišer et al., 2018).
The conclusions, however, should be considered with some care.We foresee three possible caveats that advise a caution in interpretation and application of these results.First, low number of samples per species implies sampling deficiency.In our sampling scheme we detected common species with the largest ranges, and single-site endemics.If the sampling was not sufficient, species with small ranges cannot be told apart from species with small to medium sized ranges.The species with small to mid-sized ranges might modify frequency distributions of range sizes, and if the frequency distributions of surface species changed differently from frequency distributions of subterranean species, the conclusions might change.Addressing this caveat would require additional sampling, ideally using a standardized monitoring scheme to control for the sampling effort and false negative records.
The second issue is that the sampling regions do not entirely overlap.The ranges increase with the latitude, a phenomenon known as Rappaport effect (Zagmajster et al., 2014;Hof et al., 2008).Sampling of genus Gammarus populations was on average more intense in different mountain massifs, lying in northern latitudes (fig. 1) than sampling of genus Niphargus populations.Unfortunately, the number of adequately sampled species is not high enough to appropriately model the variation of range size taking into account geographic latitude and the ecological character of species.This omission, however, is less like given that the studied region covers Downloaded from Brill.com09/17/2023 10:29:08AM via free access only approximately 10 degrees of latitude, and that the region was not strongly influenced by Quaternary glaciations.The third issue is that our analysis compares only amphipod crustaceans.The results cannot be applied to other animal groups directly, given the differences in biology among animal taxa, especially their capacities of locomotion and dispersal ecology.Any generalizations of these results need to be made with care.
We recognize that the above caveats cannot be resolved within this study, hence the results of our study need to be considered only as preliminary.We, however, find it interesting that the development of molecular methods and different quantifications of biodiversity, challenged the longstanding paradigm of subterranean endemism.Systematic research in different countries and using different taxa is needed to resolve the issue.

FIGURE 1
FIGURE 1 Distribution map of analyzed gammarid and niphargid species from Iran.Red circles show gammarid distribution records.Black triangles show previous niphargid distributions.Blue triangles show new distribution of Niphargus populations in this study.
FIGURE 2 Bayesian tree of 54 specimens as inferred from the 28S and COI gene sequences.All specimens were collected from Iran, new samples fall into two unsupported clades labelled as clade A and B. Numbers above line indicate the Bayesian posterior probabilities.Downloaded from Brill.com09/17/2023 10:29:08AM via free access

FIGURE 4
FIGURE 4 Violin plots, representing maximum linear extent for Gammarus and Niphargus species delimited using ABGD (A) and bPTP (B).The white dot in the middle indicates the median value in the range size, the thick black bar in the center represents the interquartile range and the thin black line shows 95% of all data.The distributions of surface species are not larger than the distributions of subterranean species (Mann-Whitney U test, p = 0.794 (A), p = 0.837 (B)).MOTUs found on a single localities are excluded.