Linking Hydrogeology and Ecology in Karst Landscapes: The Response of Epigean and Obligate Groundwater Copepods (Crustacea: Copepoda)

: Groundwater invertebrate communities in karst landscapes are known to vary in response to multiple environmental factors. This study aims to explore the invertebrate assemblages’ composition of an Apennine karst system in Italy mainly described by the Rio Gamberale surface stream and the Stiffe Cave. The stream sinks into the carbonate rock and predominantly feeds the saturated karst into the cave. For a minor portion, groundwater ﬂows from the epikarst and the perched aquifer within it. The spatial distribution of the species belonging to the selected target group of the Crustacea Copepoda between the surface stream and the groundwater habitats inside the cave highlighted a different response of surface-water species and obligate groundwater dwellers to the hydrogeological traits of the karst unit. Our results suggest that fast endorheic inﬁltration routes promoted the drift of epigean species from the surface to groundwater via the sinking stream while most of the obligate groundwater dwellers come from the perched aquifer in the epikarst from diffuse inﬁltration pathways.


Introduction
Groundwater invertebrate communities in karst landscapes are known to vary in response to multiple environmental factors over a wide range of spatial and temporal scales [1][2][3][4][5][6]. The link between surface water in the recharge area and karst groundwater is supposed to drive species composition in groundwater environments, despite the intrinsic difficulties in understanding the spatial distributions of water flow and storage between the surface and the underground compartments and how surface-water and groundwater species respond to the morphology, hydrology, and the speleogenesis evolution of the karst system [7].
A mature karst aquifer shows a heterogeneous spatial distribution of hydraulic conductivities in the range of 10 −10 -10 −1 ms −1 . The coexistence of low conductivities of the rock matrix (10 −10 ms −1 ) and the highest ones of large conduits draining the system (10 −1 ms −1 ) [8][9][10] are reflected in geomorphological features acting as preferential pathways of intensive groundwater circulation, with turbulent flow in karst conduits, and low water speed in the capacitive subsystems.
Karst hydrodynamic depends on the organization of the karst network, which is both highly heterogeneous and difficult to characterize. Many studies addressed the functions of the karst vadose zone and in particular its importance with respect to the epikarst [11]. The near-surface weathered zone of exposed carbonates at the rock-soil interface constitutes 2 of 18 a "recharge" zone for a karst system. Mangin [12,13] introduced the term "epikarst" to denote this zone and a perched aquifer within it at the top of the vadose zone [14]. Later, the concept was refined by Pipan and Culver [15,16] and Kozel and Pipan [17] under an integrated ecohydrological perspective due to the role of epikarst invertebrates for their potential role as tracers of water movement from the recharge area to the karst groundwaters [2,15,18].
The narrow fissures show a limited flow, are intermittently active through the hydrogeological year and constitute important zones of water and organic matter storage, working as "living chambers" for small-sized invertebrate species [1,5,6,19,20]. In some cases, a karst aquifer may be fed by either the diffuse low-flow infiltration of the epikarst and an allogenic recharge through surface runoff which drains large areas of insoluble rock or low permeability soils where surface water flows directly to adjacent soluble carbonate bedrock. In the latter case, the recharge of a karst aquifer occurs along sinking or losing stream channels via infiltration of surface water through porous streambed sediments or fast endorheic infiltration routes of streams flowing in the recharge area [21,22].
Under a hydrological perspective, recharge areas of karstic aquifers are subordinated to infiltration mechanisms which mainly determine aquifer discharge [23]. Ground watersurface water exchange is related to many factors such as geology, hydrological conditions and landscape alterations [24] and can be both spatially and temporally dynamic [25]. Open streamflow in karst very often disappears underground and emerges again in different karst formations in the subsurface, such as springs and active caves [26]. The flow regime in open karst streamflow relies mostly upon the interaction between the ground water and the surface water. They are connected hydraulically through numerous karst forms governing the water exchange between the surface and the subsurface karst formations [27,28]. Poulain et al. [11] identified two types of flow regimes in the karst. The diffuse infiltration flow is the transfer of water through rock patches with low permeability and it is identified as the slow system of recharge. Contrarily, the vertical infiltration flow represents the transfer of a high amount of water through fast transmission routes like fractures, sinkholes, sinking streams and high permeability layers. The shaft flow regime in a cave is also connected to water infiltration through the vadose zone [1,29]. These complex interactions between surface water and groundwater fluxes may affect the spatial distribution of organisms in the different cave microhabitats, the ratio between epigean species (i.e., stygoxenes which accidentally or occasionally enter groundwater through fast and slow infiltration pathways and coming from surface water bodies) and obligate groundwater dweller species (i.e., stygobites which complete the whole life cycle in groundwater). Hence, groundwater communities may be composed of epigean species coming from the sinking stream, or from the diffuse water recharge via the epikarst and the vadose zone, and of stygobite species whose evolutionary history has taken place in the small fractures of the epikarst or the saturated zone of the karst aquifer. This dual origin makes it difficult to disentangle the different routes each species may have potentially followed to reach and settle permanently or occasionally in groundwaters.
Despite the good knowledge of the groundwater dependence on surface-water discharge and physico-chemistry [18,[30][31][32][33], poor information is available on the effects of such connections on the composition of the groundwater assemblages. Whether species' drift, both epigean and obligate groundwater dwellers, at the outlets of karstic aquifers has been analysed in detail [5,6,[34][35][36][37][38][39][40][41]. However, the reverse situation, i.e., the dynamics of species' colonization from the recharge area to different groundwater habitats into a cave is poorly known, and only a few contributions have experimentally addressed this issue [4,40]. As regards the unsaturated zone of the karst, Pipan and Culver [15], Moldovan et al. [42], Liu et al. [43], and Pipan et al. [44] studied how some invertebrate species were distributed in various cave aquatic habitats in the unsaturated karst while others were found only in a single dripping pool.
In karst landscapes, it could be hypothesised that almost all the species living in the benthic layers of a sinking stream feeding an aquifer could be passively flushed out from their native habitat, thus entering groundwater in the karst, with special regard to surface-water species see [35,45,46]. It may also be argued that hyporheic species, either stygoxenes or stygobites, may be drifted during the high-discharge period of the streams flowing across the recharge area and following the endorheic flux may enter groundwater, together with water, organic matter, and eventually pollutants [47]. By the way, they may work as hydrological tracers of such connections [15].
The purpose of this study was (i) to explore species distribution patterns across a hydrological continuum of a stream flowing in the recharge area of a karst aquifer fed by the same stream; (ii) to assess which species were able to disperse from the surface to the underground; (iii) to evaluate which species may have followed an alternative route for the colonization of the saturated karst or the dripping pools in the unsaturated karst. To this end, we analysed the spatial species composition of a surface low order stream (Rio Gamberale stream) and the groundwater habitats in an active karstic cave directly fed by the stream and by the diffuse infiltration in the recharge area reaching the epikarst. The two systems are hydraulically linked by a sinkhole which allows the water of the Rio Gamberale stream to supply the Stiffe Cave. We selected the Crustacea Copepoda as the target taxon because they are the most abundant and species-rich meiofaunal group in groundwaters [48][49][50]. Copepods are also known to be good natural tracers of groundwater flowpaths [5,6,15,39,40].

Study Area
The Stiffe Cave (SC) (695 m a.s.l.) is a karst complex located on the north side of the Ocre Massif, a mountain chain 21 km far from L'Aquila (Abruzzo, Italy) [51]. The active cave is characterized by the presence of a subterranean perennial stream fed by the Rio Gamberale (RG) stream, which sinks into "Pozzo Caldaio" (1255 m a.s.l.) located in the "Altopiano delle Rocche", 2.6 km southwest from SC [52,53]. The RG is an Apennine low order stream of the Aterno-Pescara River basin. It originates from Vado di Pezza (1500 m a.s.l.) and crosses, north to south, the "Altopiano delle Rocche" for 10 km before reaching the sinkhole (Figure 1).
The recharge area of the karst aquifer is 51.08 km 2 . The endorheic water input of the RG (Q RG 360 L s −1 ) constitutes about 70% of the total water in the SC network (Q SC 500 L s −1 ). The remaining 30% of the cave water comes from diffuse infiltration paths in fractures of the carbonate rock filling the limestone aquifer [53].
The Ocre Mountains form a closed unit, hydraulically disconnected from the mountain ridges around. The water table is located about 90-100 m below the ground level of the "Altopiano delle Rocche" and it gives rise to a perched aquifer disconnected from the deeper saturated karst aquifer.
In the karst landscape of "Altopiano delle Rocche", there are several large and small fractures in the limestone from which water reaches the cave, across large chambers filled by groundwater that feed a subterranean stream flowing along the whole accessible section of the cave. From a stratigraphic perspective, the SC is developed into three superimposed distinct layers sub-horizontally oriented [52]: (i) white micrites with homogenous texture; (ii) fine and medium yellow scabrous textured calcarenites; (iii) yellow micrites interposed between thin layers of green clay ( Figure 1D). The stratigraphy directly influences the mineral composition of the streambed, which is composed of limestone pebbles constrained in clay, black limestone sand, mafic sand, medium-fine grained limestone sand, silt and clay. The cave has a total length of about 2.3 km, of which 1 km is accessible for tourism, and a vertical development of about 200 m (+186 m from the entrance). The remnant section of the cave is represented by permanently flooded large chambers which cannot be explored except by expert cave divers. The main drainage line (where the sampling survey was carried out) is characterised by seasonal floods which periodically inundate the galleries during the aquifer recharge period.

Sampling Survey
A stratified random sampling procedure was adopted for sampling the stream and the cave habitats. In RG, 21 sites were sampled: 5 from the benthic habitat (coded RG_b1-RG_b5) and 16 from the hyporheic zone (coded RG_h1-RG_h16). In SC, a total of 16 sites were sampled: 2 from dripping pools (SC_cp1 and SC_cp2), 5 from the benthic habitat (SC_b1-SC_b5), and 9 from the hyporheic habitat along the subterranean stream flowing into the cave (SC_h1-SC_h9). Except for RG_e5 and RG_h16 which underwent drought in summer, we sampled each site twice a year, in winter (in December 2014 and January 2015) and in spring/summer (in May and June 2015).
The hyporheic habitats were sampled using a Bou-Rouch pump [54]. For each sample, 10 L of water were pumped at 50 cm depth and filtered through a hand net (mesh size 60 µm). The surface benthic samples were collected with a Hess sampler (mesh size: 60 µm; depth: 10 cm) and the dripping pools inside the cave were sampled with the aid of a hand net (mesh size: 60 µm) (Figures 1-3).    For each biological sample, a set of chemical and physical parameters (temperature, pH, O2 expressed in mg L −1 , and electrical conductivity in µS cm −1 ) was measured on the field with a WTW 3430 SET G multi-parameter probe. Biological samples were preserved in 80° ethyl alcohol and taken to the laboratory where the specimens were sorted under a Leica M205C stereomicroscope and identified to species level. The Crustacea Copepoda were selected as the target taxon for the analyses because they represented more than 70% of the collected individuals of 14 invertebrate taxa. Copepods were identified to the species level and categorized in stygobites (SB) and stygoxenes (nSB). and weighted. The dry-weighted samples were burnt at 540 • C (4 h) in a muffle furnace and re-weighted to determine POM amount as the difference between dry and ash mass [55] (Supplementary File_Environmental Variables).
Biological samples were preserved in 80 • ethyl alcohol and taken to the laboratory where the specimens were sorted under a Leica M205C stereomicroscope and identified to species level. The Crustacea Copepoda were selected as the target taxon for the analyses because they represented more than 70% of the collected individuals of 14 invertebrate taxa. Copepods were identified to the species level and categorized in stygobites (SB) and stygoxenes (nSB).

Environmental Variables
Differences of environmental parameters between RG and SC sites were investigated by Principal Component Analysis (PCA) prior data standardization in z-scores. Means were used in the analysis. For the PCA only the parameters for which the standard deviation (SD) was different from zero were retained for the analysis. In addition, for the correlation coefficient between variables > 0.95, only one of the two collinear parameters was retained.

Biological Variables
Before running any type of analysis, we assessed the exhaustiveness of the sampling effort in RG and SC through species richness estimators. To assess how much larger the total number of species in both systems can get through repeated sampling we used one parametric estimator (Michaelis-Menten; MM) and five non-parametric ones (Chao1, Chao2, Jackknife 1, Jackknife 2 and Bootstrap) [56,57]. Each estimator calculates the potential species richness (S) in function of the sample size (Sobs). For each estimator, a curve of the evolution of the S predictor was obtained by gradually increasing Sobs. Values were estimated with the E-PRIMER software by means of 999 randomizations without replacement.
To assess species compositions in the benthic habitats and the hyporheic zone of Rio Gamberale stream and in the groundwater habitats inside the Stiffe Cave, for each sampling site, the sum of the species abundances of temporal replicates was calculated. After that, abundance data were converted as incidence data (presence/absence matrix) to ensure comparability of samples differing in size and sampling method. Differences between assemblages were investigated by firstly computing a distance matrix between each possible pair of sites. We selected as distance coefficient the Watson, Williams, and Lance's coefficient (WWL) [58] which is the one complement of the Sörensen index [59]. The WWL coefficient was also chosen because it is a simplification of the Bray-Curtis coefficient when it is applied to incidence data. It was calculated as follows: where j and k are the compared sites; a is the number of species co-present in j and k; b and c are respectively the numbers of exclusive species in sites j and k. Assemblage distances and sites clustering were highlighted through a Principal Coordinates Analysis (PCoA). Following the procedure suggested by Clarke [60], we tested the significance of the groups detected by PCoA using a one-way permutational analysis of variance (PERMANOVA; [61]). The significance level α was set at 0.05 (9999 permutations). A separate test for homogeneity of dispersions using the PERMDISP routine was done before performing PERMANOVA [62]. A post-hoc analysis between groups was conducted with pairwise permutation t-tests. We assessed single species contribution to the between-group overall average dissimilarities (OAD) with a similarity percentage analysis (SIMPER). The output SIMPER profile (empirical profile) was then compared with three theoretical SIMPER profiles distributions obtained with a permutation procedure (PER-SIMPER) proposed by Gibert and Escarguel [63]. The PER-SIMPER algorithm can generate distinct permutation profile distributions starting from different null hypotheses (H0). For each H0 a slightly different permutation algorithm is applied.
The first SIMPER permutation profile distribution was built under the null hypothesis (H0 1 ): the species distribution among assemblages is driven only by niche-assembly processes. This implies that the species distribution is controlled exclusively by niche numbers and breadths [64]. The second SIMPER permutation profile distribution was built starting from the null hypothesis (H0 2 ): the species distribution among assemblages is controlled only by the passive dispersal of the species [65]. A cut-off at a cumulative dissimilarity of 70% was selected.

Environmental Variables
Fourteen environmental variables with a SD = 0 (Supplementary File_Environmental Variables) were selected for the PCA and listed in Table 1. Table 1. Summary of the environmental parameters (mean ± SD) measured in the sampling sites of RG and SC (LOD = limit of detection) with SD = 0.  The PCA biplot ( Figure 4) explained 47.3% of the total variability of the environmental parameters and clearly divided the RG and the SC sites. The RG sites exhibited greater variability compared to the SC sites. The RG sites were described by higher pH values, O 2 concentration, lower Na + and Clconcentration. The RG hyporheic sites showed higher variability in SO 4 2content. Some RG hyporheic sites (RG_h12, RG_h13, RG_h14) showed higher levels of ammonium (from 2.32 to 2.81 mg L −1 ). SC sites were generally characterized by lower pH values, higher concentrations of K + and DOC and lower POM concentration. SC benthic sites were more similar to each other in their physico-chemistry than the SC hyporheic sites.

RG
The hydrochemistry of the phreatic waters into the cave mirrored, for several parameters, the one measured in the RG stream, differing in the higher DOC concentration and lower POM concentration measured in SC groundwater. The oxygen concentration was higher in the subterranean waters, and the upstream sites close to the spring feeding the RG stream. The temperature showed a tendency toward higher values in RG, even if more variable in the hydrological year, while in the cave groundwater it was lower, and less variable over time. All the other measured parameters resulted below the limit of instrumental detection.  The hydrochemistry of the phreatic waters into the cave mirrored, for several parameters, the one measured in the RG stream, differing in the higher DOC concentration and lower POM concentration measured in SC groundwater. The oxygen concentration was higher in the subterranean waters, and the upstream sites close to the spring feeding the RG stream. The temperature showed a tendency toward higher values in RG, even if more variable in the hydrological year, while in the cave groundwater it was lower, and less variable over time. All the other measured parameters resulted below the limit of instrumental detection.
The comparison between the SIMPER empirical profile and the PER-SIMPER model distribution of species ( Figure 6) highlighted that species distribution in all cases was primarily driven by niche-assembly processes, regardless of the dispersal potential of taxa. For groups A and B, A and C, and B and C the empirical SIMPER profiles were enveloped in the H0 1 PER-SIMPER distribution model for 100%, 62.5% and 96.15%, respectively.  The comparison between the SIMPER empirical profile and the PER-SIMPER model distribution of species ( Figure 6) highlighted that species distribution in all cases was primarily driven by niche-assembly processes, regardless of the dispersal potential of taxa. For groups A and B, A and C, and B and C the empirical SIMPER profiles were enveloped in the H0 1 PER-SIMPER distribution model for 100%, 62.5% and 96.15%, respectively.

Discussion
All eogenic caves have a few basic components. Water enters the subterranean system at the rock-soil interface, which typically has many small fractures and cavities with complex horizontal and vertical pathways, e.g., the epikarst [16]. Eventually, water percolating by the vadose zone reaches dripping pools in the drier section of active caves or directly feeds by means of diffuse infiltration the saturated (phreatic) zone which may be represented by perennial streams, subterranean lakes, and springs [40,45,66].
The Stiffe Cave has two different kinds of hydrological connections to the surface, the main water input is represented by the Rio Gamberale stream located in the recharge area, and a secondary contribution by diffuse infiltration water via the epikarst. Although significant differences in hydrochemistry were observed between RG and SC in the relative concentration of POM and DOC, and in the overall variation in temperature, the fastflowing of the surface stream water did not allow an enrichment in Ca 2+ , Na + and K + , suggesting that the main drains are represented by fast-flowing conduits [5,6,32]. The downstream sites of the RG are the only sites with higher concentrations of SO 4 2-, NH 4 + and Cldue to the presence of a wastewater treatment plant and intensive agriculture. In groundwater, these markers of pollution decreased likely because of the good potential of self-purification along the water infiltration pathways together with the relatively low concentrations of these pollutants in the surface flowing waters. The connectivity between overlying or adjacent surfaces and caves have received great attention since the seminal contribution of Rouch [34], pre-dated by a systematic sampling of the Baget karstic system, with the first comparative analysis of the stygoxene and stygobite harpacticoid assemblages found at the main entrance and at the outlet of the Baget karstic system [34,67]. Simões et al. [68] addressed a similar approach covering several karst caves and sampling both the phreatic zone and the epikarst in tropical Central Brazil. Moreover, Moldovan et al. [42] showed the effects of habitat fragmentation in shaping subterranean metacommunities and how habitat connectivity influenced the dispersal of crustacean fauna in an active cave of the Western Carpathians (Romania).
Despite an initial impetus in studying the faunal response to the connectivity between the recharge area and the groundwater environments in caves, the phreatic zone has received minor attention, despite the increasing studies of epikarst biodiversity, likely for the intrinsic difficulties in sampling the saturated karst in caves. Most studies of the biodiversity in the saturated karst were developed at the main outlets of such aquifers, karstic springs and boreholes having the primacy [5,20,37,39,40,45,69,70].
Our study is intended to provide further insights on the hydrologic dynamics governing the "biological connectivity" between the recharge area and the karst groundwater underneath. Copepod assemblages' compositions in Stiffe Cave were significantly different from the ones from the Rio Gamberale stream. A double signal emerged by PCoA because a clear divide was observed between the hyporheic sites of the subterranean stream (phreatic waters) and the concretional dripping pools vs. the hyporheic sites of the RG, whereas some overlapping was observed in the benthic invertebrate communities of the surface stream Rio Gamberale and the groundwater habitats of the Stiffe Cave. Group B was represented by sites from the benthic layers of both RG and SC where non-stygobite species prevailed, and the presence of stygobite species was limited to the cave sites. The non-stygobite species, shared by the benthic layers of the two groups of sites, are more easily drifted from sites of the RG stream and enter passively together with the stream water into the saturated karst of the SC. The species that have best shown this behaviour were B. pygmaeus, B. echinatus, A. crassa, and C. staphylinus. Group A, almost exclusively composed by the hyporheic sites of the SC, was best represented by the presence of the stygobite species Diacyclops paralanguidoides and Diacyclops sp. 1 which were exclusive to these sites together with Pesceus schmeili, exclusively found in SC benthic and hyporheic habitats. Moreover, higher incidence of Diacyclops paolae and Parastenocaris crenobia in hyporheic sites than in benthic sites of the SC accounted for the arrangement of this group. Group C clustered the hyporheic sites of RG stream in which the stygobite species Eucyclops intermedius and Diacyclops clandestinus were found; the former being exclusive of the hyporheic habitat of the RG stream, the latter was also found in one site only in SC. The presence of the stygobite species Diacyclops paolae, Parastenocaris crenobia, and to a lesser extent, Attheyella paranaphtalica and Elaphoidella phreatica, contributed to this grouping.
The general biological pattern observed suggests that the dissimilarity between the surface recharge area and the subterranean waters was dictated by the different composition in stygobite species, which are clear indicators of two different kinds of groundwaterdependent ecosystems: the hyporheic zone of RG and the hyporheic zone of the subterranean stream in SC. Conversely, stygoxene species represent the elements in common between these connected environments due to their higher tendency to passive dispersal, especially during the high discharge of the surface stream. In summary, the high hydrological connectivity between the surface and the underground stream is not reflected in a shared stygobite fauna.
The PER-SIMPER analysis was run to infer the ecological dynamics leading to such patterns, allowing comparison between profile distributions without focusing on single taxa or their rank in the profiles. Based on this consideration, the copepod assemblages of the hyporheic sites of the RG (Group C) and SC (Group A) showed profiles mostly described by niche-assembly, and taxa distribution is predominantly controlled by the number and breadth of niches available in each assemblage. In group C, niches for stygobites are fewer and/or narrower. In group A, a lower number of niches for surface-water species are available, thus determining a very low potential of replacement of stygobites by the occasional surface-water species which enter the saturated karst. In the benthic sites of RG (Group B) niches were more numerous and wider than in the hyporheic sites of the cave stream (Group A) and in the hyporheic sites of RG (Group C) (Supplementary file_SIMPER AvsB/SIMPER AvsC/SIMPER BvsC).
The epikarst of the SC is well developed, and dripping is consistent also in summer, but it is difficult to sample directly because the cave is active with alternation of large chambers, falls, perennial pools with more or less standing waters, and running waters, thus defining a complex saturated karst network. In the few terraced dripping pools either stygobite and stygoxene species were sampled, namely: Elaphoidella phreatica, E. plutonis, A. paranaphtalica, the stygoxene species P. schmeili, B. echinatus, M. viridis, B. pygmaeus, and B. zschokkei (Supplementary File_Incidence), and none of these species were exclusively found in this habitat.
The origin of the obligate groundwater copepods collected into the cave waters remains open to question. The species collected in the cave habitats and not found in the intensively sampled surface stream, leave room for the hypothesis that they may have colonized the subterranean waters via the epikarst. The water that fills the epikarst can infiltrate through different voids (i.e., fractures, conduits or shafts), and the route followed by the water is difficult to trace [15]. As such, the epikarst acts as a transitional zone between epigean and subterranean environments, where organic matter and other resources are transported and redirected. The species fall from above, carried by the ubiquitous dripping into the network of cavities where the environmental filters offer ideal conditions for a few species [16,31,42,44]. It is not unlikely that the "primary filter" may be dimensional, and only some stygobites, with the dominance of the worm-like and small-sized Harpacticoida Parastenocarididae and some Canthocamptidae [71], may have been incorporated in the epikarst drops, coming from the subsurface perched aquifer which works, as we say here, as a "stygobiological storage zone", and settled permanently in the saturated karst of the SC. The primarily reduced body size of the copepods [48], even smaller in groundwater species [72], makes them good candidates for dripping from the cave ceiling [15].
Our results suggest that fast endorheic infiltration routes promote the drift of epigean species from the surface to groundwater, while most of the obligate groundwater dwellers may have originated from diffuse infiltration pathways feeding a subsurface perched aquifer which fuels the epikarst.

Conclusions
If knowledge of the hydrological connectivity between surface-water and groundwater is no longer questionable, how such connections are reflected in different biological assemblages' compositions in both compartments is less known. The ecological dimension of the subterranean environments is, indeed, underestimated both for regulatory reasons and the hidden status of their biodiversity, not directly visible to most and dominated by the small-sized bacteria, protists, and invertebrates [73,74]. Further, they are considered ecologically less important than other living forms. So, why protect the invisible? The consequence of this perception is that research on the minute groundwater biodiversity is systematically neglected and less funded [75]. While much has been done to advance the knowledge of subterranean biodiversity, it remains still hindered by the difficult access to the habitats where many invertebrate species permanently live. Despite that, the groundwater biodiversity provides key ecosystem services for humans' well-being and should be adequately monitored and protected. For this reason, the spatial resolution of management actions should match as closely as possible the scale of relevant ecological patterns and processes which may affect species distributions, covering both the recharge area and the underlying aquifer, especially in the highly vulnerable karst landscapes. In fact, biological entities should be thought of in terms of dynamic processes governed by species-environment interactions, taking also into account the niche breadth of the species [76], competition, and predation among coexisting species in time and space [77,78]. On the other side, the potential for active, non-random dispersal of the species, which drives the dispersal-assembly hypothesis, may be important, but it seemed to not regulate the assemblage compositions in our case of study. In this context, new insights on species niche breadths [79][80][81] are essential for protecting groundwater biodiversity. Indeed, human activity can entail shifting in niche breadths, thus increasing the extinction risk of this unique fauna.  Acknowledgments: San Demetrio ne' Vestini municipality is acknowledged for authorizing access to the Stiffe Cave. The authors appreciate Gruppo Speleologico Aquilano and Gruppo Grotte e Forre "F. De Marchi"-CAI L'Aquila for their support in exploring and sampling the cave.

Conflicts of Interest:
The authors declare no conflict of interest.