Geochemical and microbial factors driving crustacean assemblages in adjacent aquifer units within the same aquifer

Aquifers harbor unique and highly adapted species, contributing to critical ecological processes and services. Understanding the key factors driving invertebrate assemblages in aquifers is a challenging task that, traditionally, has primarily been achieved in karst. This study aimed to uncover the factors influencing the composition and functionality of groundwater crustaceans (dimensional range from 0.036 to 1 mm) in a volcanic aquifer in central Italy. The aquifer consisted of three adjacent aquifer units (AUs) showing different geochemistry (i.e., sulfate-depleted, K-rich and, alkaline earth). We adopted a multidisciplinary approach, integrating hydrogeology, geology, microbiology, and ecology to determine whether the environmental differences that we highlighted in the three AUs were reflected in the biological assemblages. We unveiled significant differences in both the taxonomic and functional composition of groundwater crustaceans across the three AUs, and these patterns remained consistent throughout the survey period. Notably, the sulfate-depleted AU lacked groundwater-obligate species, burrowers, and stenothermal and moderately stenothermal species. The K-rich and alkaline-earth AUs had different species; however, these species exhibited similar functions related to locomotion, diet, and feeding habit. Stenothermal and moderately stenothermal crustacean species were only found in the K-rich AU, which lacked epigean species. Our findings suggest that major ions (SO${}_{\mathrm{4}}^{\mathrm{2}-}$, Ca²⁺, NO${}_{\mathrm{3}}^{-}$, and K⁺), trace elements (B, Al, V, Se, and Ba), microbial factors, and carbohydrate catabolic profiles might be the main descriptors of groundwater-obligate species abundances in the volcanic aquifer. Our findings revealed a correlation between the abundances of groundwater-obligate crustaceans and low-nucleic-acid (LNA) cells, suggesting a potential selective feeding behavior of groundwater invertebrate species on the aquatic microbial community. Our research emphasizes the need to consider diverse hydrogeological contexts within individual aquifers. Potential avenues for future research should further consider food web dynamics in groundwater communities and their impact on carbon and nutrient cycling.

1 Introduction

The Earth's largest unfrozen reserve of freshwater is found within aquifers, geological formations responsible for storing and transmitting groundwater (Ferguson et al., 2021). They harbor highly specialized microbial and metazoan species adapted to life in permanent darkness and under low-energy conditions (Malard et al., 2023a). Groundwater microbes and metazoans play crucial roles in providing ecosystem services, such as carbon recycling, pathogen removal, pollutant bioremediation, sediment mixing, and burrowing (Fillinger et al., 2023; Mermillod-Blondin et al., 2023). The variation in the taxonomic and functional composition of the invertebrate assemblages of a given aquifer results from a complex interplay of abiotic (e.g., hydrochemistry, porosity, and flow velocity) and biotic factors (e.g., species competition and cross-kingdom interactions), acting on both evolutionary and ecological scales (Zagmajster et al., 2023). However, disentangling the combined effects of these taxonomic and functional drivers remains a difficult task (Malard et al., 2023b).

Aquifers can be subdivided into specific adjacent sections, known as aquifer units (AUs), based on geological and hydrological criteria. A species' ability to spread across different AUs depends on extrinsic factors (e.g., the degree of connectivity among AUs) and intrinsic characteristics of the species, which are mainly defined by the suitability of their functional traits (e.g., body size, shape, locomotion, and feeding habit) and the ability to adjust or adapt these traits in response to the new habitat template (Fišer et al., 2023). If a species lacks trait plasticity, it is likely to struggle to establish itself in a new environment. This phenomenon is known as habitat (or environmental) filtering (Cornwell et al., 2006; Kraft et al., 2015). Groundwater crustaceans, which are the dominant metazoans in global groundwater ecosystems (Stoch and Galassi, 2010), are known to undergo environmental filtering due to abiotic factors such as pore size, flow rate, connectivity to the surface, and hydrochemistry (Trontelj et al., 2012). Environmental filtering is known to hinder the distribution of groundwater crustaceans in adjacent AUs. The organisms tend to accumulate at the boundaries of the AUs and appear unable to disperse across them (Iannella et al., 2020; Vaccarelli et al., 2023). The role of trophic filters, which relate to the availability and type of food resources and encompass both environmental filtering and competition, has been less investigated. Groundwater crustacean species in a given aquifer compete for both space and food resources (Bregović et al., 2019). Recent criticisms have pointed out challenges with respect to accurately distinguishing between environmental filtering and competitive interactions, as both can produce similar patterns in species distribution (Cadotte and Tucker, 2017). Traditionally, research on the effect of environmental filtering on groundwater crustaceans has primarily focused on karst aquifers (Maurice et al., 2016; Culver et al., 2021). However, there is a growing recognition of the importance of non-karst aquifers, such as volcanic aquifers (Wurzbacher et al., 2020). This has resulted in an intensified effort to gain a better understanding of groundwater ecosystems across a wide range of hydrogeological contexts (Hahn and Fuchs, 2009).

In this study, we aimed to unravel the factors shaping groundwater crustacean assemblages within a water table aquifer where groundwater circulates in volcanic and sedimentary formations. Preliminary analyses carried out 2 years prior to this study indicated the occurrence of adjacent, chemically different AUs within the same aquifer (Ademollo et al., 2012). We adopted a multidisciplinary approach – combining hydrogeology, geology, microbiology, and ecology – to determine whether the environmental differences that we highlighted in the three AUs were also reflected in the biological assemblages. More specifically, our goal was to pinpoint the primary factors that influence the taxonomic and functional composition of crustacean assemblages across the three AUs, including the following: (i) one hydrogeological factor (water table depth), used as proxy for the isolation of the AUs from the surface (Aquilina et al., 2023); (ii) physical–chemical factors, encompassing temperature, pH, dissolved oxygen, and major ions and trace elements; and (iii) microbial factors, such as microbial cell abundance and metabolic profiles. This study ventures into the relatively uncharted territory of volcanic aquifers and addresses critical gaps in our understanding of the main factors shaping the groundwater invertebrate assemblages.

2 Materials and methods

2.1 The hydrogeological aquifer settings and sampling survey

We conducted our study on the eastern flank of the Monti Sabatini aquifer, near the valley of the Tiber River, covering an area of 30 km² in the northeastern province of Rome (Italy). We focused on the Sabatini volcanic complex, in particular its eastern border, where the volcanic products overlap the Pleistocene sedimentary layers. The sampling site is in a rural area across the Fosso Fontanalarga, a minor right-hand tributary of the Tiber River, with small–medium villages, sparse houses, agricultural fields, and a few industrial activities, including a quarry for the exploitation of lapideous volcanic products (namely the “Via Tiberina Yellow Tuff” formation; Lombardi and Meucci, 2006). Pleistocene K-alkaline volcanic products of the Roman co-magmatic province (Sottili et al., 2010), mainly pyroclastic products (such as ignimbrites and tuffs), overlap in angular unconformity in Pliocene and Pleistocene sedimentary deposits of marine and continental origin (Fig. 1).

Figure 1 Geological map of the study area and sampling point locations. Head contours referring to the 2014 survey are plotted (heads in meters above sea level).

Holocene alluvial deposits can be found in the valley thalwegs. The water table aquifer is hosted in the volcanites, with groundwater circulating mostly in the porous and fractured lithoid levels of the ignimbrites and lavas, and in the marine and continental sands and gravels below. Low-permeability levels of the Pliocene marine clays and silts form the bottom of the aquifer and isolate it from the deep geothermal reservoir hosted in the thick sequence of the Meso-Cenozoic fractured limestones below (Parrone et al., 2020). The Pleistocene pyroclastic products, herein referred to as volcanites, typically exhibit a medium permeability. As the volcanic formations outcrop more extensively (Fig. 1), aquifer recharge occurs mainly from these rocks compared with the underlying sedimentary layers. Groundwater, mainly exploited for irrigation and to supply households not connected to the public drinking water supply system, flows at depth from northwest to southeast through both volcanic and sedimentary layers, with the latter being more important towards the edge of the volcanic system (Preziosi et al., 2013; Parrone et al., 2019).

2.2 Field and laboratory geochemical analysis

Preliminary data, collected in March 2012 (Ademollo et al., 2012), indicated the occurrence of three adjacent AUs with distinct hydrogeochemical features, characterized by sulfate-depleted, alkaline-earth, and K-rich waters, respectively. In this study, we conducted two sampling campaigns: the first was carried out between November and December 2014, while the second took place from October to December 2015. During these campaigns, we collected 20 groundwater samples from a total of 10 wells (Fig. 1; Table A1). The distribution of the wells was as follows: (i) wells QA7, QA21, QA24, and QA30, situated on the western side of the aquifer, were associated with the AU characterized by sulfate-depleted groundwaters; (ii) wells QA8, QA9, and QA29, situated along the valley of Fosso Fontanalarga, were within the AU with alkaline-earth groundwaters; (iii) wells QA13, QA15, and QA25, located on the eastern side, belonged to the AU characterized by K-rich groundwaters. For each well, we recorded essential field data, including GPS coordinates, elevation and depth of the well, depth to the water table, oxidation–reduction potential (ORP), temperature, pH, dissolved oxygen, and electrical conductivity using multiparametric portable probes (Aquaread, AP-2000). To ensure accurate measurements, we purged the wells until the physical–chemical parameters stabilized before sampling. Water samples for chemical analyses were filtered on-site with 0.45 µm membrane filters and stored in polyethylene bottles treated with 1 % HNO₃. For major cations and trace metal determination, samples were stabilized by adding HNO₃ at pH 2. Alkalinity was determined in the laboratory by HCl titration on 50 mL of sample within 24 h of sampling. Anions were determined by ion chromatography (IC; Dionex, DX-120), major cations were determined by inductively coupled plasma optical emission spectroscopy (ICP-OES; PerkinElmer, Plasma 400), and trace elements were determined by inductively coupled plasma mass spectrometry (ICP-MS; Agilent Technologies, 7500c). Samples for dissolved organic carbon (DOC) were filtered at 0.7 µm with pre-treated glass-fiber filters, and DOC levels were determined using a Shimadzu TOC-5000 A analyzer. Further analytical details are reported elsewhere (Preziosi et al., 2014; Parrone et al., 2022).

2.3 Microbial and crustacean community analyses

Samples for the microbial and crustacean community analyses were collected from all wells after purging. Aliquots (2 mL) were fixed with formaldehyde (2 % final concentration), stored at 4 °C, and analyzed within 1 week of sampling using an A50-Micro flow cytometer (Apogee Flow Systems, Hertfordshire, England), equipped with a solid-state laser set at 20 mV and tuned to an excitation wavelength of 488 nm. Total cell counts (TCCs) were quantified using side scatter and green fluorescence measurements. The intensity of green fluorescence emitted by SYBR-positive cells facilitated the distinction of cell groups based on their nucleic acid content, namely cells with a low and high nucleic acid content (LNA and HNA cells, respectively) (Amalfitano et al., 2018).

Additional aliquots (15 mL) were collected, stored in sterile Falcon tubes, and promptly incubated within a 6 h time frame to evaluate the functional diversity and metabolic preferences of the bacterial communities by measuring their degradative activity toward a set of organic carbon sources using a Biolog™ EcoPlate assay (Biolog, Inc., Hayward, California, USA). The assay consists of a 96-well microplate that contains 31 common carbon sources categorized into five compound groups (i.e., carbohydrates, polymers, carboxylic acids, amino acids, and amines), as well as a control well, in triplicate. The provided substrates are labeled with the respiration-sensitive tetrazolium dye, which is reduced to formazan when microbial degradation occurs. This is a colorimetric reaction with the typical purple color development, and optical density (OD) is detected at a wavelength of 590 nm using a multi-plate reader (PerkinElmer VICTOR™ X3). OD values were measured immediately after sample inoculation (T0) and after 24 h (T24) of incubation under dark conditions at 20 °C. Microbial metabolic fingerprinting was expressed as the percentage of each single class of compounds compared to the total OD, calculated as a fraction of the total absorbance of a single sample in the plate (Melita et al., 2019; Preziosi et al., 2019). The absorbance is proportionally related to the substrate degradation.

The sampling methodology for groundwater crustaceans was based on the manual published within the framework of the European PASCALIS (Protocols for the assessment and conservation of aquatic life in the subsurface) project, which is considered the gold standard for groundwater biomonitoring in Europe (Malard et al., 2002). Specifically, we pumped 1000 L of groundwater from each well using the submersible pumps already installed and filtered the extracted volume through a 63 µm sieve. Given the small size of most groundwater crustaceans (<1 mm in length; Malard et al., 2023a), the potential concern of the animals being minced by the pump was minimized. The 1000 L samples were taken during purging of the wells. The method is designed to capture organisms present in the well and its immediate surroundings, thereby providing a representative sample of the groundwater fauna. All boreholes were fully screened at their bottom, and the diameter of the sand filter holes was appropriate for the size of the collected animals (<1 mm). The collected samples were preserved in a 70 % ethanol solution and transported to the laboratory for further analysis. In the laboratory, we sorted the samples using a stereomicroscope at 16× magnification, identifying the taxa down to the lowest possible level of morphological taxonomic units, following established taxonomic references (e.g., Dussart and Defaye, 2002). Each specimen was categorized as either groundwater-obligate species (which complete their whole life cycle in groundwater) or epigean species (which occasionally occur in groundwater but lack adaptations for a permanent dwelling). Subsequently, we evaluated four functional traits (locomotion, diet, feeding habit, and thermal tolerance), each defined by two or more categories. The locomotion trait was defined by three categories: burrowers (organisms that dig and create burrows in sediments), interstitial (organisms that live in the spaces between sediment particles), and swimmers (organisms that are adapted to swimming in the water column). The diet trait was described by two categories: fine sediments + microorganisms (organisms that feed by ingesting fine sediment particles and the attached microorganisms) and living microphytes (organisms that feed on living algae). The feeding habit trait was described by three categories: deposit-feeders (organisms that consume organic matter settled on the substrate), collectors (organisms that gather particles from the water column), and grazers (organisms that feed on biofilms or algae growing on surfaces). Finally, the thermal tolerance trait was described by three categories: eurythermal (organisms that can tolerate a wide range of temperatures), moderately stenothermal (organisms that tolerate a moderate range of temperatures), and stenothermal (organisms that tolerate a narrow range of temperatures). Each trait was assigned based on the existing literature. We attributed each species to only one category per trait, except for the amphipod Niphargus sp., which showed a range of feeding habits and thermal tolerance (Fišer, 2019). As we did not identify the amphipod specimens down to the species level, we chose to share their percentages equally per trait category to reflect this trait variability. To obtain the trait profile of each sample, we calculated the percentage of each trait category based on its abundance within the sample (Di Lorenzo et al., 2021). This means that, for each trait (locomotion, diet, feeding habit, and thermal tolerance), we determined the proportion of individuals in each category relative to the total number of individuals in the sample.

2.4 Statistical analyses

To confirm the occurrence of the three AUs, as assumed from our preliminary investigation in 2012, we performed permutational analyses of variances (PERMANOVA; Anderson et al., 2008) based on physical–chemical factors measured in situ (water table depth, temperature, pH, dissolved oxygen, and electrical conductivity), major ions, and elements in trace, respectively. We conducted two-way PERMANOVAs, incorporating an “aquifer unit” factor with three levels (sulfate-depleted, alkaline earth, and K-rich) and a “year” factor with two levels (2014, 2015), using resemblance matrices based on Euclidean distances calculated from normalized data. To ensure the robustness of our analysis, we examined potential multicollinearity among the variables using Draftsman's plots prior to PERMANOVAs. Variables exhibiting high collinearity ($\left|r\right|\ge$ 0.95 in correlation) were considered to convey essentially identical information and were consequently removed to prevent redundancy in the analysis, in line with recommendations in Anderson et al. (2008). The variables retained for the analyses served as proxies for those that were eliminated (Anderson et al., 2008). We excluded oxygen saturation from the analyses due to > 99 % linear correlation with dissolved oxygen. In line with best practices, we applied permutation of residuals under a reduced model and employed Type III sum of squares with 999 permutations. This approach offers high statistical power and more accurate Type I error control for multi-factorial, unbalanced designs (Anderson et al., 2008). When appropriate, we conducted permutational post hoc t tests.

We evaluated the adequacy of the biological sampling effort in the whole volcanic aquifer by examining the accumulation of the total number of different observed species (Sobs) as samples were progressively added (Magurran, 2021). To assess the potential increase in species richness (S) with repeated sampling, we applied five nonparametric (Chao1, Chao2, Jackknife1, Jackknife2, and Bootstrap) and one parametric (Michaelis–Menten) estimators (Magurran, 2021). We computed the estimators at each stage as new samples were added, resulting in the generation of six curves illustrating the progression of the S with increasing sample size. We conducted the analyses through 999 randomizations without replacement (Magurran and McGill, 2011). We repeated the analyses for each AU, separately.

To investigate potential differences in the taxonomic and functional composition of the crustacean assemblage and microbial community among the three AUs, we applied the PERMANOVA design previously outlined for the environmental variables. We log(x+1)-transformed the crustacean abundances and microbial cell counts before generating Bray–Curtis resemblance matrices, while the percentage of data pertaining to metazoan functional traits and HNA and LNA cells remained untransformed (Anderson et al., 2008). We added a dummy variable equal to 1 to all data to allow the analysis of values equal to 0 (Clarke and Gorley, 2005). We chose not to analyze the abundances of each crustacean species due to the diminished interpretive accuracy when the abundances of individual taxa constituted less than 4 % of the total abundance (Clarke and Gorley, 2005). In our study, only four species met or exceeded this specified threshold (Sect. 3.2). While PERMANOVA does not inherently necessitate explicit assumptions regarding the distributions of the original variables, we opted to conduct a Levene test using the PERMDISP routine prior to all analyses. We focused on PERMANOVA outcomes that were not influenced by bias due to variance heterogeneity (Anderson et al., 2008). To provide a comprehensive overview of the significant outcomes derived from the PERMANOVAs, we utilized box plots (when considered insightful for visualization).

Finally, to assess the main hydrogeological, physical–chemical, and microbial factors that influenced the composition and functionality of crustacean assemblages across the three AUs, we employed distance-based linear models (DistLM) based on the Bray–Curtis resemblance matrix (Legendre and Anderson, 1999). We conducted both conditional and marginal tests. Conditional tests involved fitting one factor after another. We applied the BEST procedure (Legendre and Anderson, 1999) to construct models utilizing the best factor combination. We assessed the AICc (Akaike's information criterion corrected for small sample sizes; Hurvich and Tsai, 1993) value for all possible combinations of predictor variables, with the smallest AICc value indicating the most suitable model. Additionally, we used R² to evaluate the proportion of explained variation in the multivariate models. The significance of the marginal tests was determined by computing the p values through permutations rather than traditional tables (Legendre and Anderson, 1999). For each test, we employed 999 permutations to obtain p values testing the null hypothesis of no relationship, either for individual variables in isolation or in a conditional context (Legendre and Anderson, 1999). Factors that individually (marginal models) or together with others (conditional models) explained > 65 % of the variance in the taxonomic and functional structure of the crustacean assemblage were considered robust, following the criteria established by Korbel and Hose (2011, 2017) and Di Lorenzo et al. (2020). We applied this cutoff criterion to prevent unreliable and exaggerated claims about scientific phenomena (Kimmel et al., 2023). We performed distance-based redundancy analyses (dbRDA; Legendre and Anderson, 1999) to visualize the ordination of the fitted values from the most robust models.

Significance levels (α) were set at 0.05 for all permutational tests, as they provide an exact test of each individual null hypothesis of interest (Anderson et al., 2008). All analyses were performed using PRIMER 6.0 and PERMANOVA+ software (Anderson et al., 2008). Box plots were generated using the ggplot2 and gridExtra libraries and the R software (R Development Core Team, 2021).

3 Results

3.1 Environmental factors and microbial community patterns

The analyses revealed no significant differences in the water table depth among the three AUs or between the 2 years (Tables A2a and S1). However, PERMANOVA uncovered differences in temperature, pH, and electrical conductivity (Tables A2a and S1). In detail, the sulfate-depleted AU exhibited a mean temperature exceeding that of the other two AUs by about 1 °C (Table 1). The alkaline-earth AU showed the lowest pH and the highest electrical conductivity values (Table 1; Fig. 2). However, no discernible distinctions were observed between the 2 years (Table S1).

Ammonium and phosphate exhibited concentrations below the instrumental detection limit of 0.05 mg L⁻¹ in all samples and AUs and were, therefore, not included in the statistical analyses. The three AUs exhibited significant differences in the major chemical components, with variations attributed to SO${}_{\mathrm{4}}^{\mathrm{2}-}$, HCO${}_{\mathrm{3}}^{-}$, Na⁺, K⁺, and Si (Tables A2b and S1; Fig. 2). In detail, the alkaline-earth AU displayed the highest HCO${}_{\mathrm{3}}^{-}$ and Si concentrations, the K-rich AU showed the highest mean concentration of Na⁺ and K⁺, and the sulfate-depleted AU displayed the lowest mean values of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and HCO${}_{\mathrm{3}}^{-}$ (Tables 1 and S1; Fig. 2). Significant distinctions between the 2 years were only observed for Si and DOC, both of which were higher in 2014 compared with 2015 (Tables A2b and S1). Significant differences among the three AUs were also noted for trace elements, namely in the concentrations of Li, B, Rb, and U (Tables A2c and S1). The K-rich AU showed the highest concentrations of Li, B, and Rb, while the sulfate-depleted AU showed the lowest concentrations of U (Table 1 and Fig. 2). No significant differences emerged between the 2 years for trace elements.

Table 1 Physical–chemical and microbial characteristics of the groundwater samples in the three aquifer units (AUs) of the Monti Sabatini volcanic aquifer. The water table depth, field parameters, major components, trace elements, and microbial community properties are reported as mean values. Superscript letters (a, b, and c) indicate significant differences (p<0.05; permutational post hoc t tests) among AUs. An asterisk (^*) indicates statistical differences between 2014 and 2015.

Total counts of microbial cells ranged from 0.5 to 9.5×10⁴ cells mL⁻¹, with the lowest abundances occurring in the sulfate-depleted AU (Table A2d and Fig. 2). A similar pattern was observed for the percentages of HNA and LNA cells, with the lowest percentages of LNA cells found in the sulfate-depleted AU, which was correspondingly richer in HNA cells compared with the K-rich and alkaline-earth AUs (Tables 1 and A2d; Fig. 2). The microbial pattern was consistent over the 2 years (Table S1).

Figure 2 Box plots showing significant differences (p<0.05; permutational t tests) in environmental and biological parameters among the three aquifer units (alkaline earth in green, K-rich in pink, and sulfate-depleted in purple) of the Monti Sabatini aquifer. Box plots are arranged into five rows: field parameters, major chemicals, trace elements, microbes, and crustaceans. Each box plot includes median lines, 25th–75th percentile boundaries, whiskers for 10th–90th percentiles, and black dots indicating outliers. BUR denotes burrowers and EUR denotes eurythermal.

We obtained significant PERMANOVA outcomes for microbial catabolic profiles related to carbohydrate, polymer, and amino acid degradation among the three AUs, while no differences across the AUs were detected for carboxylic acid and amine utilization. However, the results for polymers and amino acids were affected by variance heterogeneity (Table S1). Additionally, catabolic profiles based on carbohydrates showed inconsistencies over the 2 years of investigation. Overall, the catabolic profile of the microbial assembly of the three AUs seemed to be mainly based on carbohydrates, polymers, and carboxylic acids, with variation over time, while the profiles based on amines and amino acids were much less represented (Table 1).

3.2 Crustacean assemblages

Out of the 20 biological samples analyzed, 8 (40 %) did not contain any invertebrate specimens (Table A2e). We identified 210 crustacean specimens (<1 mm in average length), belonging to 9 species (Table A2e). Seven of these taxa were groundwater-obligate dwellers (94 % of the total abundances; Table A2e) showing typical morphological adaptations to dark environments (Fig. 3).

Figure 3 Groundwater-obligate crustacean species. The figure showcases two groundwater-obligate crustacean species found within the Monti Sabatini aquifer: (a) the harpacticoid Parapseudoleptomesochra italica (0.73 mm in length) and (b) the syncarid Meridiobathynella sp. (0.8 mm in length). The specimens exhibit the typical morphological traits of groundwater fauna: blindness, depigmentation, and hypertrophy of sensory appendages.

Approximately 60 % of the samples contained a single species only (Tables 1 and A2e). The Copepoda class (Crustacea) dominated with seven species, including five groundwater-obligate species, six species belonging to the order Harpacticoida, and one species belonging to the order Cyclopoida. At the aquifer level, three nonparametric estimators (out of five) and the parametric one indicated that 100 % of the expected invertebrate biodiversity was collected during the sampling survey. However, the remaining estimators suggested that we likely missed a minimum of 9 % (Chao2) to a maximum of 21 % (Jacknife1) of the expected biodiversity (i.e., from one to three more species; Fig. S1). Our analyses suggested that this missed biodiversity should be found mainly in the alkaline-earth and K-rich AUs (Fig. S1).

PERMANOVA analysis indicated significant differences in the composition of crustacean assemblages; however, the outcome was affected by variance heterogeneity. On the other hand, the analyses revealed significant unbiased differences in the abundances of groundwater-obligate species among the three AUs (Table S1), with the sulfate-depleted AU lacking groundwater-obligate species (Tables 2 and S1). Notably, the K-rich AU lacked epigean species. Remarkably, there were no significant differences observed between the 2 years (Table S1). Overall, the K-rich and alkaline-earth AUs exhibited comparable abundances and species richness (Table 2).

Table 2 Taxonomic and functional composition of the crustacean assemblages in the three aquifer units (AUs). Taxonomic abundances and traits are reported as the number of individuals per AU. Superscript letters (a, b, and c) indicate significant differences among AUs (p<0.05; permutational post hoc t tests). An asterisk (^*) indicates groundwater-obligate species. For Niphargus sp. (Nsp), 0.5 individuals were counted as deposit-feeders and 0.5 individuals were counted as collectors, reflecting its dual functional feeding habit traits, while 0.5 individuals were counted as eurythermal and 0.5 individuals were counted as stenothermal, reflecting its dual thermal tolerance traits. S-D denotes sulfate-depleted and E-A denotes alkaline earth. Acronyms of taxa are reported in brackets. The numbers in the reference (Ref.) column refer to the following publications: 1 – Galassi (2001), 2 – Galassi et al. (1999); 3 – Galassi et al. (2009); 4 – Galassi and De Laurentiis (2004); 5 – Dussart and Defaye (2002); 6 – Hose and Stumpp (2019); 7 – Schminke (1974); and 8 – Fišer (2019).

However, the taxonomic composition in the two AUs was different. In detail, the groundwater-obligate harpacticoid Parapseudoleptomesochra italica was the most abundant species, accounting for 67 % of crustacean biodiversity in the aquifer (Table A2e). The groundwater-obligate P. italica and Acanthocyclops agamus were found in both alkaline-earth and K-rich AUs (Table 2). Nitocrella psammophila, Pseudectinosoma reductum, and Niphargus sp. (groundwater-obligate species) occurred in the K-rich AU, while Parastenocaris sp. and Meridiobathynella sp. (groundwater-obligate species) were collected in the alkaline-earth AU, where the epigean Elaphoidella gracilis (epigean species) was also present (Table A2e). The epigean species Moraria poppei was the only species found in the sulfate-depleted AU, with a single individual collected from well QA21 in 2014 (Table A2e).

Most of the species collected in this study were burrowers and eurythermal (Table 2). The locomotion and thermal tolerance traits exhibited variation among the three AUs, primarily driven by the higher percentages of eurythermal species in the alkaline-earth AU compared with the other two AUs (Tables S1 and A2f; Fig. 2). Notably, moderately stenothermal and stenothermal species were only collected from the K-rich AU and never from the other two AUs (Table A2f). PERMANOVA outcomes indicated significant differences in the diet trait among the AUs, although this result was biased by variance heterogeneity (Table S1). The deposit-feeding category was the most prevalent among various feeding habits, with collectors and grazers being less common, but the analyses did not identify any significant differences across the AUs (Tables 2 and S1). The pattern was consistent over time for all functional traits (Tables 2 and S1).

3.3 Linear models

3.3.1 Hydrogeological factor

The linear models based on the water table depth as a descriptor accounted for ≤ 25 % of the variance in the taxonomic and functional composition of the crustacean assemblages in the aquifer (Table S2).

3.3.2 Physical–chemical factors

Concerning field factors, the multivariate linear models accounted for < 65 % of the variance in the taxonomic and functional composition of the crustacean assemblages in the aquifer, with electrical conductivity, pH, and ORP being the main descriptors. None of the marginal models explained > 65 % of the variance (Table S2).

Multivariate models based on major ions explained < 65 % of the variance in the taxonomic composition and locomotion, diet, feeding habit, and thermal tolerance traits (Table S2). However, four major ions (SO${}_{\mathrm{4}}^{\mathrm{2}-}$, Ca²⁺, NO${}_{\mathrm{3}}^{-}$, and K⁺) together accounted for 91 % of the variance in groundwater-obligate species (Table S2), with abundances significantly decreasing in samples depleted of SO${}_{\mathrm{4}}^{\mathrm{2}}$ (Fig. 4). None of the major chemical elements individually explained more than 65 % of the variance (Table S2). Concerning trace elements, multivariate linear models explained < 65 % (Table S2) of the variance in the locomotion and thermal tolerance traits. On the contrary, they explained > 65 % of the variance in the taxonomic composition (best model 71 %: V, Cr, As, Se, and U), abundances of groundwater-obligate species (best model 89 %: B, Al, V, Se, and Ba), diet trait (best model 81 %: Li, B, Cr, Ni, As, and U) and feeding habit trait (best model 71 %: Li, B, Cr, Ni, As, and U), with abundances and trait percentages increasing with increasing concentrations of all trace elements except for Cr (Fig. 4). The contribution of individual trace elements to marginal models was always < 35 %, when significant (Table S2).

3.3.3 Microbial factors

The conditional and marginal linear models based on microbial factors (TCC, LNA cells, HNA cells, and catabolic pathways) explained < 65 % of the taxonomic and functional composition of the crustacean assemblages. However, LNA cells and the carbohydrate pathways together explained 71 % of the abundances of groundwater-obligate species (Table S2), with abundances increasing with increasing LNA cells and microbes using carbohydrates as substrates (Fig. 4). The individual contribution of LNA cells to the explained variance was 66 %.

Figure 4 Distance-based redundancy analysis (dbRDA) plots showing the linear relations between the abundances (ABB) of groundwater-obligate species and the factors selected by the BEST procedure of the DistLM. The individual panels show the (a) major ions and (b) trace element and (c) microbial factors (LNA %: percentages of low nucleic acid cells; CARB: carbohydrate catabolic profile).

4 Discussion

Despite recent advancements in the field, numerous knowledge gaps persist in identifying the factors that shape the taxonomic and functional composition of groundwater crustaceans (Mammola et al., 2021). Particularly overlooked and unexplored is the analysis of groundwater assemblages within volcanic aquifers and the role of trophic factors in shaping the distribution of groundwater-obligate fauna through environmental filtering processes (Saccò et al., 2019). To confront these challenges, our study adopted a comprehensive, multidisciplinary approach that harnessed the collective knowledge spanning hydrogeology, geology, microbiology, and ecology.

While the water table depth did not vary significantly across the aquifer, our analyses confirmed and supported the occurrence of three adjacent AUs, each characterized by distinct hydrochemical facies and microbial community patterns. The sulfate-depleted AU showed the highest mean temperature; the lowest concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, HCO${}_{\mathrm{3}}^{-}$, and U; the lowest microbial abundances and percentages of LNA cells; and the highest percentages of HNA cells. The alkaline-earth AU showed the lowest pH values and the highest electrical conductivity and HCO${}_{\mathrm{3}}^{-}$ and Si concentrations, while the microbial abundances were comparable to those in the K-rich AU. Finally, the K-rich AU presented the highest concentrations of Na⁺, K⁺, Li, B, and Rb.

We examined the crustacean assemblages of the aquifer, revealing significant variations in both the taxonomic and functional composition across the three AUs. In detail, the sulfate-depleted AU lacked groundwater-obligate species, burrowers, and stenothermal or moderately stenothermal species. The K-rich and alkaline-earth AUs, which showed comparable abundances and species richness, were characterized by different species, although these species showed the same functions concerning locomotion, diet, and feeding habit. Notably, stenothermal or moderately stenothermal crustacean species occurred in the K-rich AU only, although this AU was depleted of epigean species. We used a stringent cutoff criterion (Korbel and Hose, 2011, 2017; Di Lorenzo et al., 2020) to identify the potential predictors of these differences and concluded that water table depth was not a driving factor in shaping these assemblages. The taxonomic composition seemed to be mainly driven by trace elements, such as V, Cr, As, Se, and U. The main descriptors of the abundances of groundwater-obligate species appeared to be the major ions SO${}_{\mathrm{4}}^{\mathrm{2}-}$, Ca²⁺, NO${}_{\mathrm{3}}^{-}$, and K⁺; the trace elements B, Al, V, Se, and Ba; the microbial factors related to LNA cells; and the carbohydrate catabolic profile. Finally, the trace elements Li, B, Cr, Ni, As, and U seemed to be the main drivers of the diet and feeding habit traits. Water chemistry did not appear to exert a detrimental effect on the composition or functionality of crustacean assemblages in the target aquifer, consistent with previous studies (e.g., Di Lorenzo et al., 2020), except for Cr. This result may be attributed to the likelihood that the chemicals recognized as detrimental to groundwater organisms were either not detectable or existed at concentrations that were not harmful in our study. For instance, trace elements such as Ni, Zn, As, and Li are toxic to groundwater fauna at concentrations > 150 µg L⁻¹ (Di Lorenzo et al., 2023), although these concentrations were never measured in the study area. Ammonium (which exhibited concentrations below instrumental detection limit in our study) is lethal to groundwater-obligate species at concentrations > 12 mg L⁻¹, although it causes cellular and physiological damage at concentrations ≥ 36 µg L⁻¹ (Di Lorenzo et al., 2015). Nitrate causes no harm to groundwater fauna at concentrations < 100 mg L⁻¹ (Di Lorenzo et al., 2023), which are much higher than those observed in our study. Previous research has shown that groundwater species are sensitive to pharmaceutical compounds, pesticides, and bisphenol A (BPA; Di Lorenzo et al., 2023). In our preliminary study (Ademollo et al., 2012), chlorinated pesticides were not detected, and only traces of polycyclic aromatic hydrocarbons (PAHs) and polychlorinated biphenyls (PCBs) were found, slightly above the detection limits (0.01 ng L⁻¹). High tolerance to SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (which was positively correlated with groundwater species abundances in this study) is known for many groundwater-obligate species of marine origin, such as P. italica, N. psammophila, and P. reductum (Galassi, 2001). For instance, N. psammophila has been collected from the chemoautotrophic groundwater of the Frasassi cave system (Italy) where sulfate concentrations reach up to 199 mg L⁻¹ (Galassi et al., 2017), Parapseudoleptomesochra italica has been recorded from the Movile Cave (Romania) where groundwater has a high sulfate concentration, and P. reductum has been discovered in the in sulfidic groundwaters of Melissotrypa Cave (Greece) (Brad et al., 2021). While prior studies have reported synergistic toxic effects of pollutant mixtures on groundwater species (Di Marzio et al., 2018), our research seems to find no significant differences in chemical parameters among the various AUs that could account for varying toxicity levels. For instance, in the sulfate-depleted AU, which lacked groundwater-obligate species, most of the tested chemical elements exhibited lower concentrations compared with the alkaline-earth or K-rich AUs. Studies on the ecotoxicology of groundwater organisms have been limited by the life history traits of groundwater fauna that make them often unsuitable for laboratory experiments (Di Lorenzo et al., 2019). Our study offered preliminary insights into the potential sensitivity of groundwater crustaceans to chemicals that have hitherto remained untested in prior research, including elements like U and B. The negative correlation between the abundances of groundwater-obligate species and Cr is, in our view, not directly related to the toxicity of this element. Total Cr is present in the aquifer at concentrations ranging from 0.4 to 1.8 µg L⁻¹, which are considered harmless based on the available literature data (Di Lorenzo et al., 2023). However, this might be a statistical artifact, as lower abundances of groundwater-obligate species are found in the sulfate-depleted AU, where the average Cr concentration is slightly higher than in the other AUs. Therefore, we observe this negative correlation in Fig. 4, but we suggest that the absence of groundwater-obligate species in this AU is related to microbial factors, as we will explain later.

While the sulfate-depleted AU lacks groundwater-obligate species, the alkaline-earth and K-rich AUs differ in terms of species composition. This outcome seems to suggest that, although the AUs are hydrogeologically connected, the groundwater species collected in this study seemed not to migrate across the AUs. This finding aligns with the literature. Iannella et al. (2020) observed that European groundwater-obligate harpacticoid species (which represent 67 % of the species in our study) were unable to disperse across boundaries between two adjacent AUs. Accordingly, Vaccarelli et al. (2023) observed that the dispersal of groundwater-obligate copepod species in the eastern Lessinian Massif (Italy) is constrained by non-fractured igneous rocks, as it appears to be in our study. Only two groundwater-obligate species, P. italica and A. agamus, were found in both the K-rich and alkaline-earth AUs. The remaining species collected in this study were unique to one AU only, such as the groundwater-obligate P. reductum (a Tertiary relict of ancient marine origin; Galassi et al., 1999), which was collected from the K-rich AU, along with N. psammophila. The genus Nitocrella, which has direct marine origin, serves as an indicator of ancient evolutionary events (Galassi et al., 2009). Notably, the K-rich AU showed the highest concentrations of Li, B, and Rb, which are characteristic elements of the deep geothermal facies and possibly of fossil marine waters of Neogene age (Duchi et al., 2003). The presence of stenothermal and moderately stenothermal species in this AU (and only in this one) suggests that the habitat of this K-rich AU may be conservative and, therefore, suitable for preserving ancient evolutionary lineages with no close relatives in surface environments, such as P. reductum.

Two out of the six indicators revealed that a small percentage of taxonomic diversity (ranging from one to a maximum of three species) was not captured in this study. This finding is commonly encountered because groundwater sampling is essentially a blind process (Mammola et al., 2021). Wells serve as windows through which we gain insight into the subterranean biodiversity in the portions of the aquifer surrounding them, but groundwater habitats extend more extensively (Ficetola et al., 2019). However, we speculate that the absence of groundwater-obligate species in the sulfate-depleted AU might be partly due to the low percentages of LNA cells found in this AU, rather than, or in addition to, incomplete sampling. The LNA and HNA cell counts exhibited variations among the three AUs. Specifically, the sulfate-depleted AU displayed the lowest LNA cell percentages. As detected by flow cytometry, LNA and HNA cell groups are considered constitutive traits of aquatic microbial communities, typically comprising cells of varying sizes, genome content, and phylogenetic affiliations (Gasol et al., 1999; Proctor et al., 2018). Notably, HNA cells were recognized as an active fraction of the bacterioplankton community, and their abundance was reported to positively correlate with heterotrophic production rates in freshwaters of different origin (Rubbens et al., 2019). In contrast, LNA cells were traditionally thought to represent a more dormant or quiescent portion of the aquatic microbial community (Lebaron et al., 2002). They were reported as small-sized microorganisms with slower metabolic activity and a wide range of survival strategies suitable for thriving in adverse conditions, including oligotrophy (Hu et al., 2022). In previous groundwater studies, the percentage of LNA cells and the use of carbohydrates were linked at either low or high nutrient levels (Melita et al., 2019). This is explained by the role played by carbohydrates as an important energy-rich carbon source for aquatic microbial metabolism (Arnosti et al., 2014). These observations seem to suggest that groundwater quality can directly affect the functional properties of the aquatic microbial communities with implications for the pattern of the energy fluxes among organic matter, microbes, and the organisms located in the upper levels of the food web. In this study, the high linear correlation (r=0.80) between groundwater crustacean abundance and LNA cells (supported by an explained variance R²>65 %) seems to suggest that groundwater crustaceans might derive significant benefits from the presence of LNA cells. Considering the high abundances of deposit-feeders and crustaceans feeding on sediments and microorganisms in the alkaline-earth and K-rich AUs compared with the sulfate-depleted AU, we venture to speculate that they may selectively feed on LNA cells. If confirmed by further studies, our findings would suggest that the feeding behavior of groundwater crustaceans might play a role in structuring the groundwater microbial community and biomass, with potential consequences for subterranean carbon turnover and nutrient cycling.

Total cell counts also varied across the AUs with the sulfate-depleted AU showing the lowest values. However, we observed only a moderate linear correlation (r=0.5, R² = 25 %; Table S2) between the TCC and groundwater crustacean abundance, indicating that the TCC explains a limited portion of the variation in crustacean abundance. This result might seem counterintuitive, as most of the collected species are known to feed on fine sediments and microorganisms. Nevertheless, there should be a logical explanation for this outcome. The planktonic microbial community represents a fraction of the total aquifer microorganisms found in the interstitial water volume due to detachment from sediments (Flemming and Wuertz, 2019). As the crustacean species examined in this study are not filter-feeders, it is probable that they consume a portion of microbial planktonic cells smaller than the portion of the sediment-attached ones. This feeding habit could be the underlying cause of the weak correlation observed in this study. This observation raises the possibility that a more robust correlation between diet or feeding habit and LNA cells might have emerged had we included sediment-attached cells in our analysis.

On this matter, we did not observe a significant correlation between DOC and groundwater crustacean assemblages. This suggests that DOC is likely not directly utilized by the species collected in this study and rather needs to be processed by bacteria before becoming accessible to metazoans (Foulquier et al., 2010; Griebler and Lueders, 2009; Segawa et al., 2015). Furthermore, DOC concentrations did not significantly differ across the three AUs, with consistent values of < 1 mg L⁻¹, in line with the values normally observed in groundwater systems (Foulquier et al., 2010). DOC entering subterranean environments from the surface primarily comprises stable and recalcitrant components that resist bacterial degradation, leaving only a small fraction available for microbial communities (Shen et al., 2015). Our analyses indicated that the sulfate-depleted AU had the highest proportions of microbes metabolizing polymers, which are more complex and resistant to degradation than simple carbohydrates or carboxylic acids (Oest et al., 2018; Melita et al., 2019). Additionally, we observed variations in the proportions of microbes utilizing different substrates from year to year, likely influenced by surface-produced organic matter types (Saccò et al., 2019; Melita et al., 2023). Our findings prompt a discussion about the minimum quantity of energy resources needed to sustain a resident groundwater-dependent crustacean assemblage. To provide a rough estimate, considering that a prokaryotic cell contains approximately 25 fg of carbon (Griebler et al., 2002), the sulfate-depleted AU in our study had an average microbial biomass of 2.65–4 mg C L⁻¹. As no groundwater-dwelling species were collected in this AU during our survey, further studies are necessary to assess if this biomass level could be a limiting factor for groundwater-dependent crustacean assemblages. Finally, we acknowledge two potential limitations of our study: the incomplete sampling of crustacean biodiversity and the timing of the sampling survey. As we missed a small percentage of the expected biodiversity, our findings should be considered preliminary. Complete biodiversity sampling in groundwater is challenging due to the many impediments associated with the essentially blind nature of groundwater sampling (Mammola et al., 2021). Furthermore, while eDNA and metagenomics approaches currently face challenges in certain contexts, particularly with markers such as 18S (Korbel et al., 2017), classic DNA barcoding using cytochrome oxidase I gene (COI) markers has been shown to be an effective tool for delineating cryptic groundwater species (Altermatt et al., 2023). In addition, while seasonality in groundwater recharge may influence geochemical processes (Jasechko et al., 2014), previous studies have shown that seasonality has relatively little impact on the distribution of biota in groundwater ecosystems. For instance, Korbel and Hose (2015) demonstrated that habitat structure, water quality, and site attributes are the key environmental variables influencing groundwater metazoan distribution, with only minimal variance explained by seasonality. Extending the survey timing, sampling effort, and implementing classic DNA barcoding and eDNA approaches in future studies would be helpful to confirm our results and address these potential limitations.

5 Conclusions

Our multidisciplinary study delved into the taxonomic and functional composition of groundwater crustaceans and attempted to unravel the intricate dynamics of the crustacean communities in three AUs within the same aquifer. We assessed significant variations in the crustacean distribution and functional traits across the AUs, with the sulfate-depleted AU standing out as a seemingly inhospitable environment for groundwater-obligate species. Our findings pointed to the crucial role of microbial communities in driving the composition of groundwater-obligate crustacean assemblages. Additionally, the study stimulates a discussion on the sensitivity of groundwater-obligate species to aquifer settings. Our research underscores the importance of singling out diverse hydrogeological contexts within individual aquifers. Potential avenues for future research encompass metagenomic studies on specific microbial taxa that are at the base of groundwater food webs, while stable isotope analyses would help elucidate the dietary preferences and food web dynamics of groundwater-obligate crustaceans and their impact on nutrient cycling.

Appendix A

Table A1 Characteristics of the sampling sites within the three aquifer units (AUs).

^* The unit “m a.s.l.” denotes meters above sea level.

Table A2 Physical–chemical, microbial, and crustacean data (taxonomic and functional) of the groundwater samples in the three aquifer units (AUs) of the Monti Sabatini aquifer (S-D: sulfate-depleted; E-A: alkaline earth; K: K-rich). Panel (a) presents the field data (WT: water table depth, in meters above sea level; ORP: oxidation–reduction potential, in mV; T: temperature, in °C; DO: dissolved oxygen, in mg L⁻¹; EC: electrical conductivity, in µS cm⁻¹). Panel (b) provides the major components (in mg). Panel (c) shows the trace elements (in µg). Panel (d) presents the microbial community characteristics. Total cell count (TCC, in cell mL⁻¹) and low-nucleic-acid (LNA) and high-nucleic-acid (HNA) cells are expressed as percentages of total absorbance. The other abbreviations used in this panel are as follows: CARB – carbohydrates; POL – polymers; CARB_A – carboxylic acids; AM – amino acids; and AMIN – amines. Panel (e) gives the taxonomic composition information. The abbreviations used in this panel are as follows: Pit – Parapseudoleptomesochra italica; Psp – Parastenocaris sp.; Nps – Nitocrella psammophila; Pre – Pseudectinosoma reductum; Egr – Elaphoidella gracilis; Aag – Acanthocyclops agamus; Mpo – Moraria poppei; Msp – Meridiobathynella sp.; Nsp – Niphargus sp.; and SB – cumulative abundances of groundwater-obligate crustacean species. In each sample (ID), abundances are reported as the number of individuals per 1000 L of groundwater. Panel (f) presents the functional composition (as a percentage). The abbreviations used in this panel are as follows: BUR – burrowers; INT – interstitial; SWI – swimmers; FS-M – fine sediments + microorganisms; LM – living microphytes; D-F – deposit-feeders; COL – collectors; GRA – grazers; EUR – eurythermal; MST – moderately stenothermal; and STE – stenothermal.