Sampling campaigns were carried out in different periods. The Presciano spring was sampled in January 2004 (5 sites), March 2004 (5 sites), July 2004 (5 sites), September 2004 (5 sites), January 2012 (5 sites), March 2012 (5 sites), July 2012 (7 sites), and October 2014 (7 sites). The Capo Pescara-Santa Liberata (hereafter CP_SL) spring was sampled in August 2010 (6 sites), September 2010 (3 sites), March 2011 (7 sites), July 2011 (3 sites), November 2011 (7 sites), May 2012 (6 sites), September 2012 (7 sites), January 2013 (11 sites), and March 2015 (11 sites). The Cavuto spring was sampled in October 2006 (6 sites), December 2006 (6 sites), February 2007 (4 sites), March 2007 (5 sites), October 2014 (8 sites), and March 2016 (8 sites). Finally, the Chiarino spring was sampled in July 2014 (6 sites), August 2014 (10 sites), May 2015 (10 sites).

All the microhabitats detectable were sampled (karst fractures, sediments devoid of vegetation, sediments with macrophytes, submerged mosses) and, where available, running waters from the rheocrenic sites and lentic waters from the limnocrenic sites. The surveys concerned both inbenthic (i.e., the topmost sediments at a maximum depth of 10 cm) and subsurface (at a depth of 40 cm) sites for each spring. Inbenthic samples were collected with a Hess sampler (mesh size: 60 μm; diameter: 40 cm) by washing the sediments contained in the cylinder and moving the sediment particles under the current generated. The finest grains and the fauna dislodged were thus collected into the downstream net (Dole-Olivier et al., 2014). Subsurface samples were collected by hammering mobile steel piezometers with a 5 mm-hole-screened tip within the sediment at a depth of about 40 cm; they were then connected to a membrane pump according to Malard et al. (2004). Ten liters of interstitial water were extracted, filtered with a 60-μm-mesh net, and fixed in a 70% alcohol solution. After sorting the invertebrates collected, the copepods were identified to the species level and classified as stygobites or non-stygobites according to their degree of adaption to the groundwater life (e.g., miniaturization, depigmentation, anophthalmy, mature females with larger but fewer eggs than their epigean counterparts, according to Galassi et al., 2009a).

We used a multiparametric probe (WTW Multi 3430 SET G) to measure, in the field, electrical conductivity at 25°C (µS cm^-1), pH, dissolved oxygen (mg L^-1), and temperature (°C) of water collected from each sampling site at each time and spring. To exclude the presence of a chemical contaminant affecting the patterns of species richness and beta diversity, we measured, for a subset of 48 sites encompassing all the four springs, the concentrations of 108 contaminants related to anthropogenic pressures (uncontrolled municipal waste landfills, intensive agriculture and farming, urban settlements, industries) occurring in the recharge areas of the springs. The targeted chemicals were sulphates, phosphates (method: APAT CNR IRSA 4020 Man 29 2003), nitrates, nitrites, ionized ammonium (method: 135 APAT CNR IRSA 4020 Man 29 2003), metals (method: EPA 3005A 1992+ EPA 6010C/APAT CNR IRSA 3150C Man 29), pesticides (method: EPA 3510C 1996+ EPA 8270D 2007), volatile organic compounds (method: EPA 5030C 2003+ EPA 8260C 2006; EPA 8260 C 2006), and hydrocarbons (method: EPA 3510C 1996+ EPA 8270D 2007). Their concentration was measured in the laboratory by retaining a 2-L sample for each sampling site.

Total, stygobite, and non-stygobite copepod species richness were calculated for all the sampling sites and sampling dates. To assess the potential influence of in situ physical-chemical parameters (water temperature, electrical conductivity, pH and dissolved oxygen) on the temporal dynamics of species richness within the target springs, we used random intercept Generalized Linear Mixed Models (GLMMs), setting spring identity as a random effect, through the “lme4” package (Bates et al., 2015). In the mixed models, the fixed effects included the explanatory variables (and possibly their interactions) presumed to directly affect the response variable, while the random effects accounted for the portion of data variability deriving from inter-site relatedness along unmeasured dimensions (e.g., geographic proximity, Zuur et al., 2009). Different sites from a single spring were likely to be more similar in terms of biotic assemblages than sites from different springs; therefore, setting spring identity as random effect allowed us to model environment-richness and environment-beta diversity (see below) relationships taking this autocorrelation structure into account.

Before fitting the models, we applied the data exploration protocol according to Zuur et al. (2010). The distributions of the explanatory (i.e., physical-chemical parameters) and response (i.e., richness values) variables were inspected, and possible outliers in the dataset were detected and removed when they were likely to derive from data recording errors. Then, a stepwise Variance Inflation Factor (VIF) analysis was performed through the “usdm” package (Naimi et al., 2014) to lower multicollinearity among the explanatory variables, which could lead to incorrect estimation of the coefficients of the fixed effects and related p-values (Zuur et al., 2010; Cerasoli et al., 2021). Specifically, VIF = 3 was set as the threshold VIF score for variables’ retention (Zuur et al., 2010).

Once the set of non-collinear physical-chemical parameters was identified, an initial Poisson GLMM was fitted on these parameters, also including the corresponding interaction terms. To facilitate model convergence, the values of the physical-chemical parameters were scaled (i.e., s c a l e d   X i = X i − μ σ ) before model fitting. Then, model selection was performed by sequentially dropping, in turn, one of the fixed effect terms (starting from interactions) and testing the statistical significance (considering p-value ≤0.05 as threshold) of the withheld term through a χ 2 -based likelihood ratio test, as well as by looking at the reduction in the corrected Akaike’s Information Criterion (AICc) score of the corresponding model (Zuur et al., 2009).

Successively, we performed a simulation-based model validation on the refined GLMM (i.e., the model obtained at the end of the variables’ selection phase). Specifically, we used the “simulateResiduals” function from the “DHARMa” (Hartig, 2022) package to simulate new data from the fitted model and successively derive scaled residuals by comparing the simulated data with observed values. This process was repeated 500 times and the resulting Q-Q plot and scatterplot of scaled residuals versus model predictions were examined to confirm uniformity of residuals and the absence of overdispersion. In case of remaining uniformity or overdispersion issues, we fitted again the GLMM on the same set of predictors by using a negative binomial distribution instead of Poisson distribution. The simulation-based validation was then applied to the newly obtained model to verify the absence of any remaining issues.

After the model selection and validation steps, we used the “sjPlot” package (Lüdecke, 2022) to extract from the final model: (i) estimated value, standard error, 95th percentile confidence interval and associated p-value for each fixed effect term; (ii) variance attributed to the random effect; (iii) Nakagawa’s pseudo-R² (Nakagawa et al., 2017); (iv) AICc and log-Likelihood scores; (v) marginal effect plots for the fixed effect terms resulting as statistically significant (i.e., p-value ≤ 0.05).

We used the “betapart” (Baselga et al., 2022) package to investigate within-spring temporal variation in taxonomic and phylogenetic beta diversity. Total beta diversity and its two components, turnover and nestedness-derived diversity (hereafter nestedness), were first measured using the corresponding multiple-site dissimilarity metrics as in Baselga (2012). Calculation of the beta diversity components was performed for each spring and sampling date through the “beta.multi” function and using the Sørensen family of dissimilarity indices (i.e., BetaSIM for turnover, BetaSNE for nestedness, and BetaSOR for total beta diversity).

To assess the statistical significance (considering α = 0.05 as threshold for Type I error) of the observed variation in beta diversity over time, we implemented a random resampling approach to derive a distribution of values for each of the above-mentioned diversity metrics, computed on random subsets of the available sites. This analysis was carried out, for each spring, using the “beta.sample” function. Specifically, 30 iterations of the resampling procedure were run, considering only the sampling dates for which at least 7 different sites were sampled: in this way, 35 unique combinations of four distinct sites were extracted for each of the considered sampling dates; 30 out of these 35 unique combinations of sites were then randomly chosen to compute multiple-site beta diversity metrics. Statistically significant differences among the distributions of dissimilarity values obtained for the different sampling dates were assessed, for each diversity metric, using non-parametric tests because preliminary linear models did not fulfill the normality, homoscedasticity, and independence assumptions (Zuur et al., 2009). For the Cavuto, Chiarino, and Presciano springs, only two sampling dates included enough sites to perform the resampling procedure. Therefore, for these springs we carried out two-sided Mann-Whitney U tests to compare the distributions of resampled beta diversity values from the two considered dates. For the CP_SL spring we had five sampling dates (March 2011, November 2011, September 2012, January 2013, and March 2015) with at least 7 sampled sites. Thus, we first performed a Kruskal-Wallis Rank Sum Test for each diversity metric, to check for statistically significant differences among the five sampling dates. Then, to compare each pair of sampling dates, we performed post hoc two-sided Mann-Whitney U tests with Benjamini and Hochberg (1995) p-value correction.

The phylogenetic counterparts of the above-mentioned beta diversity metrics were computed through the “phylo.beta.multi” function of the “betapart” package. As a proxy of phylogenetic information to be included in this calculation, we derived an ultrametric phylogram of the copepod species recorded in the studied springs, based on Linnean taxonomy (Supplementary Figure S1), through the “ape” (Paradis and Schliep, 2019) package. The tree was rooted in the Harpacticoida order, following previous studies indicating the Harpacticoida as more basal than the Cyclopoida (Eyun, 2017; Bernot et al., 2021).

To assess if, and to what extent, the physical-chemical conditions of the springs influenced variation in the components of taxonomic and phylogenetic beta diversity, we initially fitted GLMMs with beta distribution (as beta diversity metrics are bounded between 0 and 1) and logit link functions through the “glmmTMB” (Brooks et al., 2017) package. The candidate explanatory variables included were the maximum, minimum and standard deviation values of temperature, pH and dissolved oxygen computed, for each spring and sampling date, across the sampled sites. This choice permitted contrast of the beta diversity observed at each sampling date with the corresponding inter-site variability and extreme values of the physical-chemical parameters. Within the resulting set of “summary” physical-chemical parameters, the ones not affected by multicollinearity were selected through the same VIF-based approach used for the species richness modelling.

As for the GLMM built for species richness, explanatory variables were scaled before model fitting. A different model was fitted for each possible combination of two explanatory variables and corresponding two-way interaction term. The sample size was relatively low (n = 26, corresponding to the spring × sampling date combinations), so that including more than two explanatory variables would have likely hampered model convergence. Nonetheless, preliminary model fitting and validation runs still showed convergence issues and unreliable estimations of model parameters (e.g., negative residual and random effect variances). Thus, we moved from the initial GLMMs to Linear Mixed Effect (LME) models with normal distribution, fitted through the “lme4” package. Among the LMEs fitted on the distinct combinations of explanatory variables, we selected the ones leading to the lowest AICc score. We performed model refinement through sequential dropping of the retained explanatory variables. Then, the simulation-based model validation as well as the extraction of model statistics and related marginal effect plots were performed following the same steps described for the analysis of species richness.

A total of 48 species, comprising 22 stygobites and 26 non-stygobites, were recorded within the 168 sites sampled across the four springs (Supplementary Table S1). Nine species were collected from at least one site in all springs, while twenty-five species were found in a single spring (Supplementary Table S2). Among these latter, it is worth mentioning some still undescribed species, such as a Nitocrella species sampled from a single pool in the Cavuto spring and two new species belonging to the canthocamptid genera Moraria and Paramorariopsis which were found in the Chiarino spring. Bryocamptus echinatus was the most widespread species in all the springs analyzed, with Bryocamptus zschokkei and Diacyclops paolae being among the top-five species in terms of occurrences in three out of the four springs (Supplementary Table S3).

Total species richness varied among the four springs, as well as within each spring over time (Figure 3). The Presciano spring showed the highest species richness (SR), with maximum values in July 2012 and October 2014 (22 and 24 copepod species, respectively). The CP_SL spring showed lower maximum SR (15 species) than the Cavuto (18 species) and Chiarino (20 species) springs. The Presciano spring harbored the highest number of stygobites (hereafter Sb) among the four springs, with 9 Sb species found in July 2012 and October 2014, followed by CP_SL (6 Sb species in March 2015), Cavuto (5 Sb species in October 2014) and Chiarino (6 Sb species in May 2015).

Total beta diversity was higher than 0.6 for almost all the sampling dates, with species turnover generally higher than nestedness (Figure 4). An apparent increase in overall beta diversity over time, driven by rising turnover, emerged for the Chiarino spring and the Presciano spring (starting from March 2012 for the latter), while the values of the three diversity components more markedly fluctuated over time in the Cavuto and CP_SL springs. For Capo Pescara-Santa Liberata spring, statistically significant differences among the five compared sampling dates emerged for all the dissimilarity metrics (nestedness: χ2 = 99.44, p < 0.001; turnover: χ2 = 94.53, p < 0.001; total beta diversity: χ2 = 52.99, p < 0.001). Fluctuations in turnover, nestedness and total beta diversity values between the compared pairs of sampling dates were statistically significant for all the springs, with few exceptions (Supplementary Tables S4, S5).

When accounting for among-sites phylogenetic differences, trends across springs and sampling dates were similar to those previously observed for taxonomic diversity (Figure 5). Indeed, phylogenetic turnover was higher than phylogenetic nestedness for most samples, thus suggesting that several sites host unique representatives of higher-order taxonomic groups, such as species belonging to subfamilies and/or genera not recorded in other sites from the same spring.

All the springs showed pristine conditions, as the average concentrations were below the LOD (Limit of Detection) for almost all the considered contaminants (Supplementary Table S6). Physical-chemical parameters varied considerably among the springs, with Chiarino and Cavuto showing lower temperature and electrical conductivity, and higher pH and dissolved oxygen concentration than CP_SL and Presciano (Supplementary Figure S3).

The electrical conductivity (VIF score = 5.49) determined most of the multicollinearity among the four physical-chemical variables. Once discarded, pairwise collinearity among the remaining variables was lower than 0.7 and VIF scores were lower than 3.

The Poisson GLMM relating total SR to the retained physical-chemical variables showed both overdispersion (DHARMa dispersion = 1.63, p-value = 0.046) and non-uniformity of residuals (one-sample two-sided Kolmogorov–Smirnov test: D = 0.15, p-value <0.001). The negative binomial GLMM fitted on the same fixed effect terms was not affected by these issues (DHARMa dispersion = 0.98208, p-value = 0.424; one-sample two-sided Kolmogorov–Smirnov test: D = 0.075, p-value = 0.328). The single physical-chemical parameters did not significantly contribute to the model, whereas the corresponding two-way interaction terms did (Table 2). Specifically, above average values for dissolved oxygen led to steeply increasing predicted SR when temperature (Figure 6A) and, more prominently, pH (Figure 6B) showed below average values (i.e., along negative semiaxis of the scaled gradient). Overall, the selected model did not explain much of the observed temporal variation in species richness; indeed, the marginal pseudo-R² was as low as 0.06, although it noticeably increased when considering the spring-related random effect (i.e., conditional pseudo-R² = 0.23).

The Linear Mixed Effect (LME) model fitted with total beta diversity (BetaSOR) as response variable showed non-negligible explanatory power even without considering the spring-related random effect (marginal pseudo-R² = 0.26). In particular, among-sites standard deviation in dissolved oxygen (hereafter O₂ sd) was positively related to predicted beta diversity (Table 3; Supplementary Figure S4).

The interaction between O₂ sd and maximum within-spring dissolved oxygen values (hereafter O₂ max) emerged as a statistically significant term in the LMEs fitted with turnover (BetaSIM) and nestedness (BetaSNE) as response variables (Table 4). Predictions from these LMEs showed a good fit to observed values of nestedness and turnover, particularly at the extremes (either low or high) of the O₂ max gradient, which is also confirmed by the relatively high marginal pseudo-R² (Table 4). In particular, when maximum O₂ values within the springs were below average (i.e., negative scaled O₂ max), high inter-site O₂ variability (i.e., positive scaled O₂ sd) produced an increase in predicted turnover (Figure 7A) with a corresponding decrease in predicted nestedness (Figure 7B); the opposite occurred, with steeper variations, when O₂ variability was low (i.e., when O₂ sd was below average).

None of the LMEs fitted on the taxonomic beta diversity components showed residual overdispersion or uniformity issues (p-value >0.05 in all the validation tests).

Outcomes of the LMEs fitted for the different components of phylogenetic beta diversity resembled those obtained for taxonomic diversity in terms of significant fixed effect terms and associated estimated coefficients (Supplementary Tables S7, S8). Particularly, the predicted positive effect of O₂ sd on total phylogenetic diversity (Supplementary Figure S5) as well as the synergistic influence of O₂ sd and O₂ max on phylogenetic turnover and nestedness (Supplementary Figure S6) recalled the patterns emerging from the LMEs fitted on taxonomic beta diversity.