Trapped in the web of water: Groundwater‐fed springs are island‐like ecosystems for the meiofauna

Abstract We investigated whether the equilibrium theory of island biogeography (ETIB) can be applied to the meiofauna of groundwater‐fed springs. We tested whether copepod species richness was related with spring area, discharge, and elevation. Additionally, five hypotheses are tested based on species distribution patterns, dispersal ability, and life‐history characteristics of several guilds (stygobiotic, nonstygobiotic, cold stenotherm, and noncold stenotherm species). Thirty springs in the central Apennines (Italy) were considered. A multimodel selection procedure was applied to select best‐fit models using both ordinary least‐squares regressions and autoregressive models. Mantel tests were used to investigate the impact of spatial autocorrelation in determining interspring similarity (ßsor), pure turnover (ßsim), intersite nestedness (ßnest = ßsor − ßsim), and matrix nestedness (measured using NODF and other metrics). Explicit consideration of spatial correlations reduced the importance of predictors of overall species richness, noncold stenotherm species (both negatively affected by elevation), cold stenotherm species, and nonstygobiotic species, but increased the importance of area for the stygobiotic species. We detected nested patterns in all cases, except for the stygobites. Interspring distances were positively correlated with ßsor and ßnest (but not with ßsim) for the entire data set and for nonstygobiotic, cold stenotherm, and noncold stenotherm species. In the case of stygobites, interspring geographical distances were marginally correlated with ßsor and no correlation was found for ßsim and ßnest. We found support for ETIB predictions about species richness, which was positively influenced by area and negatively by elevation (which expresses the size of source of immigrants). Low turnover and high nestedness are consistent with an equilibrium scenario mainly regulated by immigration and extinction. Stygobites, which include many distributional and evolutionary relicts, have a low capability to disperse through the aquifers and tend to be mainly confined to the springs where they drifted out and were trapped by springbed sediments.


| INTRODUCTION
The equilibrium theory of island biogeography (ETIB) proposed by Wilson (1963, 1967) to explain variation in species number on islands turned into one of the most productive research programs in ecology (Lomolino, 2000;Lomolino, Riddle, Whittaker, & Brown, 2010;Warren et al., 2015). Providing a simple mechanistic model based on extinction and colonization rates, the ETIB obtained an immense success even far beyond its original scope, being applied to virtually any kind of isolates (Walter, 2004). The widespread implications of the ETIB were already clear to MacArthur and Wilson themselves, who stated that the principles of island biogeography could be applied to any kind of isolated habitats, such "streams, caves, gallery forest, tide pools, taiga as it breaks up in tundra, and tundra as it breaks up in taiga" (MacArthur & Wilson, 1967: 3-4).
Despite the ETIB has reached the status of a scientific paradigm, few studies have explicitly tested its foundations (Whittaker & Fernández-Palacios, 2007), because of obvious difficulties in measuring colonization and extinction rates. Thus, most research that evoked the ETIB has been based on the assumed relationships between geographical parameters and extinction/colonization rates. For example, it is usually hypothesized that immigration rates should decrease with increasing distance of islands from the mainland and/or other islands and increase with island size (the so-called target effect, because larger islands are larger targets) (Lomolino et al., 2010;Warren et al., 2015;Whittaker & Fernández-Palacios, 2007). At the same time, it is assumed that extinction rates should decrease with increasing island size, because larger islands host larger, and hence more stable populations (Lomolino et al., 2010;Warren et al., 2015;Whittaker & Fernández-Palacios, 2007), and that extinction rates should also decrease with island proximity due to rescue effects (Brown & Kodric-Brown, 1977).
Springs are currently classified as groundwater-dependent ecosystems (Eamus & Froend, 2006), yet they are usually approached under an epigean perspective, and subsurface spring habitats have been almost completely neglected Stoch et al., 2016). This appears surprising, because groundwater-dependent ecosystems are becoming key research topics in freshwater biology and ecology (Howard & Merrifield, 2010). Springs may be considered as natural windows to the subterranean world, thus representing unique environments to investigate the ecology and biogeography of groundwater organisms which reach the springs from the aquifers where they spend the whole life cycle . However, most of the researches on spring ecology have dealt with plants and animals living in the spring head above the bottom substratum, with only a few studies investigating the spring subsurface meiofauna .
Moreover, these studies dealt mainly with the influence of physicochemical parameters on spring community structure , whereas no study has been addressed to the biogeography of the groundwater meiofauna inhabiting the springs. There is increasing evidence that springs are geographical isolates of conservation concern for the high levels of cryptic diversity of both surface-water and groundwater macroinvertebrates (Murphy et al., 2015), even if the explanatory factors shaping differentiation of cryptic species and microendemicity are still matter of debate (Juan et al., 2010;Rader et al., 2012). In this study, we aimed at testing-for the first time-the principles of ETIB to groundwater invertebrate meiofauna using explicitly stated hypotheses.
We considered two levels of diversity: species richness recorded in each spring and β-diversity (spatial turnover) of copepods (Crustacea, Copepoda, Figure 1). According to the principles of ETIB, we postulate that species richness should increase with spring area and decrease with isolation. Also, we expect a distance decay effect, that is, that interspring similarity in species composition should decrease with increasing interspring distances, which should lead to a correlation between β-diversity and geographical distances (Nekola & White, 1999).
Dynamic models, such as the ETIB, suppose that distribution patterns on islands are mainly influenced by present conditions, such as island distance to the mainland, interisland distances, and habitat heterogeneity, more than by historical factors. However, some species may display distributional patterns that are the result of past events, and relict models postulate that species distribution on islands was F I G U R E 1 The stygobiotic copepod crustacean Pseudectinosoma reductum Galassi & De Laurentiis, 1997. It was discovered in the Presciano spring system (Gran Sasso aquifer) and after recorded also from the Cavuto spring (Montagna Grande aquifer). This species is a "living fossil" of ancient marine origin, with a relict distribution and a very isolated phylogenetic position determined by past interisland and/or island-mainland connections followed by vicariance events, instead of present conditions, such as island distance to the mainland (see Whittaker & Fernández-Palacios, 2007).
The relative influence of current and historical factors varies according to the island systems and the organisms. In general, species with low dispersal capabilities show distributional patterns more strongly constrained by regional historical factors (such as geological events, paleogeographical setting, and paleoclimatical changes), whereas present local and regional ecological factors may be more important for taxa with higher dispersal abilities.
According to the ETIB, the species-area relationship is a consequence of the fact that extinction rates decrease and immigration rates increase with island area. By contrast, according to relict models, a positive species-area relationship results from area-dependent extinction or relaxation processes. Stygobites (i.e., organisms that complete their whole life cycle in groundwater) have lower dispersal capabilities in comparison with nonstygobiotic species (Galassi, Huys, & Reid, 2009;Stoch & Galassi, 2010;Trontelj et al., 2009), and their distributional patterns are recognized to result mainly from past geological events (but see Eme et al., 2015 for the role of spatial heterogeneity and productive energy).
On the basis of these considerations, we formulated the following predictions: 1. Species richness should correlate positively with spring area for both nonstygobites and stygobites. However, we hypothesize that, if stygobiotic species form relict assemblages, their richness should not be significantly influenced by other variables, such as discharge and elevation, which may be important for species which conform to the ETIB.

2.
If stygobites form relict populations, they should have idiosyncratic distributions that are not strongly influenced by current geography; thus, we expect that their species richness values should not be strongly autocorrelated. By contrast, if nonstygobites follow an equilibrium model, with species exchanges facilitated by interspring proximity, we expect a strong spatial autocorrelation in species richness values for these animals.

3.
Because of their idiosyncratic distributions, stygobiotic species should present low nestedness. By contrast, if nonstygobites are subject to continuous processes of immigration/extinction, we expect significant nestedness.

4.
Noncold stenotherm and nonstygobiotic species are known to have higher dispersal capabilities compared to stygobites, and hence, we hypothesize that they can move easily through the surface hydrogeological network at low altitude. So, we expect that noncold stenotherm and nonstygobiotic species richness are negatively related to spring altitude. Moreover, these species should be more widespread and show lower spatial turnover.

5.
Cold stenotherm species need a stable thermal regime. As higher groundwater discharge guarantees cold and stable thermal regime, we expect that cold stenotherm species correlate positively with discharge.

| METHODS
The studied springs are located in the central Apennines (Abruzzo region). Climate is continental, with average annual precipitation ranging from 700 to 1,000 mm. Springs are distributed in three main Thirty springs, fed exclusively by groundwater, were sampled across the three hydrogeological units analyzed ( Figure 2). All the investigated springs are pristine, perennial, and located in undisturbed areas. A stratified sampling has been performed in each spring in order to cover most of the environmental heterogeneity. Sampling was carried out by means of drift nets placed at the major openings of the bedrock in the case of karstic springs; alternatively, when sediments were present, the interstitial habitat was sampled by pumping 10 L of interstitial water with the Bou-Rouch method (Bou & Rouch, 1967) and then filtered with a hand net (60-μm-mesh size). Among the invertebrates collected, the microcrustacean copepods were identified to species level and used as target group being by far the most abundant and species-rich taxon among meiofaunal organisms. Values of species richness reported in Table 1 should be considered virtually complete (Galassi, Fiasca, & Del Tosto, 2011). Primary data are given in Appendix S1 in Supporting Information.

| Statistical analysis
Both the ETIB and relict models predict that species richness should increase with island area (species-area relationship) and should decrease with isolation. We recorded spring area as the spring surface permanently covered by water. For larger, typically limnocrene springs Measuring spring isolation is not a trivial procedure. Teittinen and Soininen (2015) used the proportion of terrestrial area in the catchment from the entire catchment area as a measure of spring isolation. This kind of measure does not seem appropriate for groundwater organisms. For true islands, isolation is usually measured as island distance to the mainland, which is assumed as the species source (Weigelt & Kreft, 2013). In the case of groundwater-fed springs (i.e., springs fed exclusively by groundwater, in contrast to springs fed by rainfall, snowmelt, and glacier melt), the species source (i.e., the mainland) is the aquifer. Thus, a perfect analogous for springs should be the distance to the aquifer network, which is however impossible to measure, especially for the meiofauna. The animals that compose the meiofauna live not only in the small and large fractures filled by groundwater, but, in the case of karstic aquifers, also in the annex capacitive subsystem , avoiding the large karstic conduits where water velocity is very high. Such patchiness in the distribution of the meiofauna, together with the intrinsic complexity of the aquifer network, makes it practically impossible to measure the groundwater flowpath from the aquifer to the spring. Thus, the subterranean pathways of water interconnections that can be used by the meiofauna remain unknown.
For these reasons, we opted for an indirect measure that might be viewed as a functional analogous of "distance to mainland." We considered spring discharge as an inverse "functional" proxy of islandmainland distance. This is because a higher discharge, depending on the flow rate and the length of the flowpath from the recharge area, indicates a larger volume of water that is transferred from the aquifer to the spring, so that it should increase immigration rates of copepods. Mainland size is usually not considered in island biogeography (but see Weigelt & Kreft, 2013 for an isolation metric that accounts for the coastline shape of large landmasses by including only the area of the part that extends into a certain buffer). However, in the case of springs, it is important to consider that the proportion of the aquifer that represents the source pool of organisms is not the same even for springs that are fed by the same aquifer. Specifically, while highaltitude springs receive water from a small portion of the aquifer (because of the proximity of the surface recharge area of the aquifer to the spring outlet), and thus may be likely colonized by individuals coming from a reduced source pool, basal springs are fed by a larger part of the aquifer and hence receive immigrants from a likely larger source pool.
In general, elevation is considered one the most important ecological factors shaping species richness patterns of copepod assemblages in springs, especially because stygobites are rare in high-altitude mountain springs, whereas they represent an important fraction of the copepod species richness in springs located at lower altitudes (Stoch, 2007).
For these reasons, we introduced both discharge and elevation as possible factors influencing immigration rates, and hence species richness. We used few environmental predictors in our models even if other variables, such as spring temperature, shading, conductivity, F I G U R E 2 Map of the study area (Abruzzo mountains) with indication of sampled springs. Springs codes as in Table 1. Top right inset shows location of the study area in Italy pH, concentration of chemical parameters, or springbed sediment texture, might be important in determining species richness (Galassi et al., 2011;Reiss & Chifflard, 2015). However, we preferred to select only ecogeographical variables that paralleled those used in ETIB formulation, leaving the study of other correlates to future research.
Environmental variables were always log-transformed to linearize their relationships with species richness values. The use of log-transformed area values is also consistent with the Gleason exponential model of the species-area relationship: where S is the number of species, A is the area, and c and z are fitted parameters (see Triantis, Guilhaumon, & Whittaker, 2012). For comparative purposes, we also applied the Arrhenius power model (log S = log c + z log A; see Triantis et al., 2012), which however provided a slightly poorer fit.
Area was log (x + 1)-transformed because of the presence of zero values for springs so small that their surface was not possible to be measured with sufficient accuracy. We used a multimodel selection procedure to identify best-fit models on the basis of the corrected Akaike information criterion (AICc). To investigate the impact of spatial autocorrelation on the relationships between environmental variables and species richness, we used both ordinary least-squares regression (OLS) and spatial autoregressive (SR) models (Beale, Lennon, Yearsley, Brewer, & Elston, 2010;Fattorini & Ulrich, 2012;Liechstein, Simons, Shriner, & Franzreb, 2002) as implemented in SAM (Spatial Analyses in Macroecology) v. 4.0 software (Rangel, Diniz-Filho, & Bini, 2010). Autocorrelation was evaluated using Moran's I index with 1,000 permutation to calculate p-values. We conducted analyses for total species richness and for stygobiotic, nonstygobiotic, cold stenotherm, and noncold stenotherm species, separately.
T A B L E 1 Hydrological characteristics, geographical location, and meiofaunal species richness of 30 groundwater-fed springs studied in Central Italy Lat, latitude (decimal degrees); Lon, longitude (decimal degrees); E, elevation (m); D, discharge (Ls −1 ); A, area (m 2 ); T, total number of species; Sb, number of stygobiotic species; Cs, number of cold stenotherm species. Codes refer to Figure 2.
To investigate correlation between interspring variation in species composition and interspring distances, we measured pairwise minimum surface distances between springs using QGIS (QGIS 2.12 Development Team, 2012). These distances express only the geographical proximity of the springs and not necessarily their degree of interconnection, which is influenced by the subterranean hydrology.
In fact, springs that are geographically very far might be connected by subterranean networks, and geographically close springs may be hydrologically rather isolated. However, we used surface geographical distance as a measure of isolation under the assumption that interspring closeness should facilitate the exchange of individuals belonging to species with high dispersal power, although the reverse is not necessarily true (i.e., species with low dispersal capabilities might migrate from a spring to another if they are connected by the subterranean hydrology). We used the approach of Baselga, Jiménez-Valverde, and Niccolini (2007) and Baselga (2010Baselga ( , 2012 for partitioning the overall ß-diversity (ßsor) among springs into true species replacement or pure turnover (ßsim) and nestedness (ßnest) components. In this respect, nestedness quantified the part of compositional change caused by ordered species loss, whereas pure turnover was related to the exchange in faunal composition caused by the local trade-off between species extinction and immigration. To assess whether variations in ßsor, ßsim, and ßnest were geographically structured, we correlated their values with interspring geographical distances using Mantel tests (Pearson correlation coefficient, with 10,000 random permutations to assess significance). Lack of significant correlations indicates that the considered species are unable to cross the surface landscape, whereas significant correlations may arise by both surface and underground dispersal.
To quantify the total degree of nestedness in each matrix, we used standard nestedness derived from the overlap and decreasing fill Although Z-scores can be used to assess whether a matrix is significantly nested or not, they cannot provide information about the "magnitude" of nestedness (Strona, Stefani, Galli, & Fattorini, 2013). For this, we calculated also the "relative nestedness" (RN) proposed by Bascompte, Jordano, Melián, and Olesen (2003), which is  (2000). Nestedness analyses were conducted using the NeD program (Strona et al., 2014) with 100 null matrices.

| RESULTS
We found significant spatial autocorrelations (Moran's I at the first  Table 1 and Figure 2 (a) (b) (c) (d) T A B L E 2 Parameter values, standard errors and associated probability levels of ordinary least-squares (OLS, a, b, c) and spatial autoregressive (SR, d, e, f) best-fit models for total, stygobiotic and cold stenotherm species richness, and log-transformed environmental variables Stygobiotic species richness was not spatially autocorrelated (I = .142, P I=0 = .658), which supports our Prediction #2.

Variable
For the complete data set, the OLS best-fit model included area and elevation as significant predictors (Table 2). Area, however, was not significant when the autoregressive model was applied (Table 2).
By contrast, area was the only predictor in the best-fit OLS model for stygobiotic species and was significant also in the autoregressive model (Table 2), which supports a relict model for stygobites (Prediction #1).
For the cold stenotherm species, the OLS best-fit model included area and discharge (Prediction #5), but these variables were not significant in the autoregressive model (Table 2). For the nonstygobiotic species, the OLS best-fit model included only elevation as a significant variable (Prediction #4), which however was not significant in the autoregressive model (Table 3). The OLS best-fit model for the noncold stenotherm species included area and elevation as significant variables, but only elevation was significant in the autoregressive model (Table 3).
In general, compared to the independent errors model, explicit consideration of spatial correlations reduced the importance of predictors of overall species richness, cold stenotherm species richness, noncold stenotherm species richness, and nonstygobiotic species richness, but increased the importance of area for the stygobiotic species.
Using the entire data set, both overall species dissimilarity among All nestedness indices indicated a significant nestedness for the complete data set and for the cold stenotherm species, which showed a very similar nestedness degree (Table 4). Similarly, both nonstygobiotic and noncold stenotherm species showed significant nestedness values (Table 4). By contrast, all indices but matrix temperature showed a nonnested pattern for the stygobiotic species (Table 4), in accordance with Prediction #3.
T A B L E 3 Parameter values, standard errors and associated probability levels of ordinary least-squares (OLS, a, b) and spatial autoregressive (SR, c, d) best-fit models for nonstygobiotic and noncold stenotherm species richness, and log-transformed environmental variables  F I G U R E 4 Relations between interspring distance and overall similarity in species composition (ßsor, a, d, g, j, m) and its nestedness (ßnest, b, e, h, k, n) and turnover (ßsim, c, f, i, l, o) components for total species richness (a-c) and stygobiotic (d-f), nonstygobiotic (g-i), cold stenotherm (j-l), and noncold stenotherm (m-o) species, separately

| DISCUSSION
In general, we found support for ETIB predictions about overall species richness, with both area and elevation being important predictors.
Area was included in all best-fit models, although its importance varied according to the species ecological categories that were modeled, and was generally diminished when the spatial component was considered. However, for stygobites, area was the only variable retained in the best-fit models and its importance further increased when the spatial component was included in the analysis (more than 80% of variance in the OLS model and more than 95% of variance in the autoregressive model [predictors + space]). Thus, our first prediction that stygobiotic richness should be influenced by area, but not necessarily by other variables, was largely confirmed. One of the explanations for the species-area relationship is the so-called habitat diversity hypothesis (Hart & Horwitz, 1991;Kohn & Walsh, 1994). According to this hypothesis, larger areas host more species because they have a greater variety of biotopes, larger environmental gradients and more diverse and abundant resources, etc. (Hortal, Triantis, Meiri, Thébault, & Sfenthourakis, 2009;Hortal et al., 2013;Fattorini et al., 2015).
Most of the springs considered in this study are basal karstic springs, where groundwater upwells through large and small fractures of the carbonate bedrock and in the alluvial sediments overlying the aquifer, thus offering several microhabitats for species survival beneath the springbed Stoch et al., 2016). However, the socalled passive sampling hypothesis (Connor & McCoy, 1979) can also be evoked. According to this hypothesis, larger islands have more species simply because they "sample" more individuals, and hence more species, from the total pool of immigrants (Connor & McCoy, 1979;Ricklefs & Lovette, 1999;Fattorini et al., 2015). These hypotheses are not mutually exclusive and can be evoked from both an ETIB and a relict perspective. However, in the case of stygobites, the lack of spatial autocorrelation and relationships with other variables suggests that these animals form relict communities (Prediction #2), for which the species-area relationship is not determined by current immigration/ extinction processes. Namely, we hypothesize the presence of a "filter effect" exerted by alluvial deposits on species with low dispersal capabilities and that determined a more selective immigration, compared to that of nonstygobites, which have higher dispersal capabilities. This "filter effect" on species washed out by groundwater into the alluvial deposits of the springs led stygobites to be trapped in the sediment matrix (see Galassi et al., 2014). Under this scenario, springs with larger areas may have sampled and retained more species during the "washing" process, whereas rheocrene springs with small areas, where the groundwater rapidly flows downstream, have a much lower "retention potential" for stygobites, resulting in a negligible "filter effect".
Interestingly, we found that our data were best fit by the Gleason exponential model of the species-area relationship, instead of the Arrhenius power function. It has been hypothesized that the exponential model is more suited when accumulation of new species is relatively slow within an area due to the probability of strong similarity of environmental conditions and species composition of neighboring areas (He & Legendre, 1996). As our system is composed of springs with similar environmental characteristics and species composition, with a relatively small variation in species richness (3-23 species), our results seem to support this hypothesis.
Elevation exerted a negative influence on the total species richness, on nonstygobiotic species richness and on noncold stenotherm species richness. This suggests that nonstygobiotic species and noncold stenotherm species richness values are both higher in loweraltitude springs, where these generalist species may reach the spring habitats not only from the "mainland" represented by the aquifer but also, and likely predominantly, by means of stepping-stone dispersal through the surface hydrographic network (Pringle, 2003), in accordance with our Prediction #4. Moreover the intrinsic mosaic feature of the large basal springs offers higher niche availability, thus supporting the co-occurrence of all the ecological categories analyzed in the present study. Cold stenotherm species richness was positively influenced by area and discharge, thus suggesting that immigration rates of species associated with cold water is enhanced by higher volumes of water, as expected according to Prediction #5. Higher groundwater discharge guarantees cold and stable thermal regime and may favor the development of aquatic vegetation that represents a trophic and spatial microhabitat where some cold stenotherm species live, as observed in Alpine ponds (Ilg & Oertli, 2014). Moreover, most high-altitude springs are rheocrenic habitats, where the groundwater that feeds the springs flushes out and permanence of water bodies is mainly determined by the springscape (Reiss & Chifflard, 2015) described by slope, current velocity, and groundwater discharge (von Fumetti, Nagel, Scheifhacken, & Baltes, 2006;Van Der Kamp, 1995).
Comparisons between OLS regressions and spatial autoregressive models are highly informative. We observed that, with the notable exception of area for stygobites, the use of spatial autoregressive models reduces the importance of predictors selected by OLS regressions. Of course, many environmental parameters are geographically structured, and removing the spatial component should indirectly lead to reduce their influence. Thus, the fact that elevation remained an important predictor of total species richness even in the autoregressive models indicates that this variable exerts an influence that cannot be subsumed by variations in other possible factors associated with spatial position. By contrast, area remained important for stygobites, but not for cold stenotherm and noncold stenotherm species. This suggests that the influence of area on these categories may be due to the influence of variables associated with the positional effect, such as the hydrological unit that can serve as a source pool, the effect of different histories that different parts of the study area may have experienced, or the presence of additional spatially structured environmental conditions that were not considered in this study.
In general, interspring geographical distances influenced both the overall species dissimilarity and the nestedness component, but not the pure turnover component (Prediction #3). This suggests that interspring proximity influences interspring species similarity mainly via variation in species richness; when the effect of species richness is removed, that is, when the pure spatial turnover is considered, the effect of geographical distance disappears. Thus, low turnover and high nestedness values, also returned by NODF and other metrics, are consistent with an equilibrium scenario mainly regulated by immigration and extinction. However, in the case of stygobiotic species, as originally predicted, we found only a marginally significant correlation between interspring geographical distances and overall species dissimilarity, whereas no correlation was found for the pure turnover component or the nestedness component. Moreover, the species per spring matrix of stygobiotic species was not significantly nested for all used metrics, except for T, which is however known to be prone to type I error (Ulrich et al., 2009). The lack of nestedness results from unexplained species presences or absences and thus indicates that the distribution of stygobiotic species is largely idiosyncratic and was influenced by biogeographic events different from the immigration/extinction processes affecting the other species. For example, geological and past hydrological processes may have led to alternate phases of isolation and connection of the aquifers (Bauzà-Ribot, Jaume, Fornós, Juan, & Pons, 2011), which may explain the co-occurrence or disjunct distribution of stygobiotic species in aquifers that are now isolated.
Distributional and evolutionary relicts are common among stygobites inhabiting subsurface spring environments, where these species can found conditions (such as a constant thermal regime, an oligotrophic status, and permanent darkness) that mirror those typical of their original habitat, represented by the aquifer feeding the springs where they are trapped (Galassi et al., 2009;Holsinger, 1988;Humphreys, 2000).

| CONCLUSIONS
Springs are known to serve as ecological refugia (variable in time and space) for surface-water and groundwater invertebrates (Cantonati et al., 2012), and also as evolutionary refugia (Botosaneanu, 1995;Davis et al., 2013) especially for stygobites whose primary habitat is represented by the aquifer itself hosting narrow endemic and paleoendemic taxa, the latter with disjunct distribution across springs fed by different aquifers. Whereas nonstygobiotic species show distributional patterns mainly influenced by real-time ecological constraints and which are consistent with ETIB predictions, stygobiotic species show spot distributions among spring habitats, being the only survivors of ancient lineages that become extinct in surface-water bodies, or, alternatively, as a result of fragmentation of a wider original distribution. In both cases, springs represented for stygobites conservative environments that "trapped" species originally depending on the stable groundwater environments of the aquifers below the spring cup (Galassi et al., 2009;Gibert & Deharveng, 2002).