Elucidating stygofaunal trophic web interactions via isotopic ecology

Summary

Subterranean ecosystems host highly adapted aquatic invertebrate biota which play a key role in sustaining groundwater ecological functioning and hydrological dynamics. However, functional biodiversity studies in groundwater environments, the main source of unfrozen freshwater on Earth, are scarce, probably due to the cryptic nature of the systems. To address this, we investigate groundwater trophic ecology via stable isotope analysis, employing δ¹³C and δ¹⁵N in bulk tissues, and amino acids. Specimens were collected from a shallow calcrete aquifer in the arid Yilgarn region of Western Australia: a well-known hot-spot for stygofaunal biodiversity. Sampling campaigns were carried out during dry (low rainfall: LR) and the wet (high rainfall: HR) periods. δ¹³C values indicate that most of the stygofauna shifted towards more ¹³C-depleted carbon sources under HR, suggesting a preference for fresher organic matter. Conversion of δ¹⁵N values in glutamic acid and phenylalanine to a trophic index showed broadly stable trophic levels with organisms clustering as low-level secondary consumers. However, mixing models indicate that HR conditions trigger changes in dietary preferences, with increasing predation of amphipods by beetle larvae. Overall, stygofauna showed a tendency towards opportunistic and omnivorous habits—typical of an ecologically tolerant community—shaped by bottom-up controls linked with changes in carbon flows. This study provides baseline biochemical and ecological data for stygofaunal trophic interactions in calcretes. Further studies on the carbon inputs and taxa-specific physiology will help refine the interpretation of the energy flows shaping biodiversity in groundwaters. This will aid understanding of groundwater ecosystem functioning and allow modelling of the impact of future climate change factors such as aridification.

Keywords:

Copepods – Trophic interactions – Predation – Beetles – Larvae – Amino acid analysis – Stable isotopes – Food web structure

Introduction

During recent decades, investigations of trophic webs have become a cornerstone for the interpretation of functional biodiversity in freshwater ecosystems. Within both lentic and lotic environments, macroinvertebrate food web dynamics play a key role in shaping process-level aquatic ecosystem attributes [1]. Aquatic faunal trophic characterization is usually conducted by employing the morpho-behavioural based concept of functional feeding groups (FFGs) [2]. Since its inception, FFGs have been extensively used in ecological assessments and biomonitoring studies, and have allowed detailed assessment of ecological patterns in both natural and disturbed environments [3,4,5].

However, despite the hydraulic and ecological continuum in groundwater dependent ecosystems, the subsurface ecosystem and the study of its food chain interactions have suffered from a conceptual disconnection from surficial aquatic habitats. The main reasons are attributable to methodological limitations [6,7], scarce aquifer accessibility [8] and the lack of interdisciplinary approaches [9]. Moreover, compared to surface freshwater ecosystems, groundwaters are subjected to relatively extreme environmental conditions: sparse organic inputs, lack of light and primary production, and truncated trophic webs [10,11,12,13]. Altogether, these unique conditions shape obligate subterranean aquatic communities (stygofauna) dominated by plastic and opportunistic trophic behaviours [14,15], whose categorization via feeding modes such as FFGs is constantly at risk of misinterpretation. As a result, our knowledge about how food web interactions shape groundwater ecological functioning and community patterns is fragmented [16].

Stygofauna—when present—play a key role in regulating both ecological and hydrological dynamics in aquifers [17,18]: they actively bioturbate the sediment, facilitate nutrient recycling and, in combination with microbial communities, degrade/retain contaminants. In groundwaters, carbon inputs (allochthonous dissolved organic carbon (DOC) and chemoautotrophic production) are mediated by microbes which are then grazed by basal stygofauna [19]. Organic matter (OM) is transferred along the trophic chain via prey-predator interactions. Therefore, OM inputs, microbial communities, and the stygofaunal trophic web, all shape the energy flows sustaining the subterranean biodiversity [20].

The incorporation of biogeochemical approaches (i.e. stable isotopes composition, fatty acids content, radiocarbon analysis) has recently led to re-evaluation of the archetype of poorly structured–and generalist-dominated–trophic dynamics in groundwaters [21]. These designs are leading a vital transition from purely descriptive to functionally-based investigations, providing wider perspectives to the field [22].

Carbon (δ¹³C) and nitrogen (δ¹⁵N) stable isotope analysis (SIA) is a well-established approach enabling quantitative investigation of food webs [23,24]. Since its initial application in groundwater trophic ecology, several studies have benefited from the insights provided by the study of naturally-occurring stable isotopes [25,26]. However, δ¹³C and δ¹⁵N SIA investigations on bulk material are limited by the mixing of tissues and different biochemical pathways [27]. These limitations can be addressed by the complementary or alternative use of compound-specific approaches.

δ¹³C and δ¹⁵N Compound Specific Isotope Analysis (CSIA) on amino acids allows detailed characterization of food web interactions [28], by focusing on compounds created by definable biosynthetic pathways. Single amino acids can be divided into essential (EAA) and non-essential (NEAA). Whilst primary producers (plants, algae and bacteria) biosynthetise de novo EAA from a bulk carbon pool, animals lack these enzymatic pathways and acquire EAA from their diet [29]. As a result, tracking of EAA allows carbon fingerprinting of food sources down to the base of food webs [30]. Concurrently, δ¹⁵N CSIA can distinguish between compounds reflecting the source isotopic signal, and that enriched with each trophic step, thus providing crucial information on prey-predator interactions [31]. The application of CSIA in amino acids has allowed a much more thorough understanding of food web dynamics in freshwater [32], marine [33] and terrestrial environments [34], but despite the greater potential than bulk analysis [35], this technique has yet to be applied to food web studies of groundwater environments.

This study is, to our best knowledge, the first based on the combination of carbon and nitrogen CSIA in groundwaters, and focuses on a calcrete stygofaunal community under two contrasting environmental conditions: low rainfall (LR, dry season) and high rainfall (HR, wet season). We hypothesise that different environmental conditions trigger species-specific adaptations that are ultimately responsible for distinct food web interactions. The specific objectives of this work are: 1) unravel OM incorporation trends across the stygofaunal community, 2) decipher the trophic habits of the species and elucidate prey-predator interactions, and 3) provide biochemically-based knowledge about trophic web interactions in arid zone calcrete aquifers.

Methodology

Study area and field work

The field work was carried out at a calcrete aquifer (28°41‘S 120° 58‘E) located on Sturt Meadows pastoral station, Western Australia, ~42 km from the settlement of Leonora (833 km northeast of Perth, Fig 1A).

The Sturt Meadows calcrete hosts a very shallow aquifer, located two to four metres below the surface, and accessible through bores characterised by water depths ranging from a few centimetres to ten metres. The bore grid was initially drilled for mineral exploration and comprises 115 bore holes of between 5–11 m in depth forming a 1.4 km by 2.5 km (3.5 km²) area (Fig 1A). The bores are unlined, except for about the upper 0.5 m which are lined with 10 cm diameter PVC pipe for stability, and capped [36]. Three sampling campaigns–two of them corresponding to low rainfall periods (LR) and one during the wet season (high rainfall, HR) [37]–were carried out in July and November 2017, and March 2018. More details about the sampling design, monitoring of water depth and hydrogeological background at Sturt Meadows can be found in Saccò et al. [38].

The high morphologically (Fig 1B) and taxonomically diverse stygofaunal community at Sturt Meadows comprises three sister species of subterranean beetles (Paroster macrosturtensis (Watts & Humphreys 2006), Paroster mesosturtensis (Watts & Humphreys 2006) and Paroster microsturtensis (Watts & Humphreys 2006) and respective larvae)), three species of amphipods (Scutachiltonia axfordi (King, 2012), Yilgarniella sturtensis (King, 2012) and Stygochiltonia bradfordae (King, 2012)), aquatic worms (family Tubificidae (Vejdovský, 1884)) and water mites (order Oribatida; Dugès, 1834). Within the stygobiotic meiofaunal community, two species of harpacticoids (Novanitocrella cf. aboriginesi (Karanovic, 2004), Schizopera cf. austindownsi (Karanovic, 2004) and four species of cyclopoids: Halicyclops kieferi (Karanovic, 2004), Halicyclops cf. ambiguous (Kiefer, 1967), Schizopera slenderfurca (Karanovic & Cooper, 2012) and Fierscyclops fiersi (De Laurentiis et al., 2001)) can be found.

Adult and larval stygofaunal specimens were collected by hauling a small weighted plankton net (mesh 100 μm, [36]) five times from the bottom through the water column of 30 boreholes (Fig 1A) selected by simple random sampling [38]. Stygofaunal abundance data across the boreholes are reported in S1 Table.

All biological samples were kept frozen (–20°C) in darkness until further processing in the laboratory where individual organisms were counted and identified (and consequently separated) to the lowest taxonomic level via optical microscopy and reference to specific taxonomic keys. Roots and sediment samples from the bottom of the aquifer were obtained through the stygofaunal haul netting procedure, and were separated by using sterile glass pipettes during the sorting in the laboratory according to the sampling campaign (LR or HR). Sediment samples were soaked in acid (0.1 N HCl) to remove inorganic carbon and dried at 60°C for 24 hours. Given the delicacy of the hydrological dynamics in shallow calcretes [39], extensive water extractions spread along the bores were avoided and preliminary tests were carried out to quantify the potential risk of dewatering the calcrete.

Bores D13 and W4 host groundwater systems which are representative of the geological conformations of the area—phreatic and vadose calcretes interspersed with clay material—and were finally selected because of their hydrological and biotic stability (lowest risk of drying and representative ranges of Sturt Meadows stygofaunal diversity) [38]. Water samples for POC (particulate organic carbon) analysis were collected using a submersible centrifugal pump (GEOSub 12V Purging Pump) after wells were purged of three well-volumes and stabilisation of in-field parameters was observed. POC samples were obtained by filtering water from the bores D13 and W4 through GF/F filters (pre-combusted for 12 hours at 450°C), washed with 1.2 N HCl to remove any inorganic carbon, and subsequently dried at 60°C for 24 hours. The field site was accessed and samples were collected with permit approval (permit number 08-003150-1) from the Department of Parks and Wildlife of Western Australia.

Sample preparation and study design

All individuals from a single taxon were pooled for each sampling campaign (LR1, LR2 or HR) and subsequently washed with MilliQ water to remove external contaminants. Subsequently, samples were oven dried at 60°C overnight and crushed to a fine powder which was stored at –20°C until further analysis (Fig 2).

Due to sample size constraints, the samples for each taxon from the two low rainfall periods were further combined. Previous metagenomics investigations, together with mesocosm experiments and field observations at Sturt Meadows provided some information about the trophic habits of beetles and amphipods [40]. Adult subterranean beetles had active predatory feeding on epigean amphipods and copepods (including group feeder behaviours) together with scavenger habits (and potential active predatory pressures) on sister species. Beetle larvae (third and last instar) showed opportunistic predatory habits with a range of prey from copepods and amphipods to adult beetles from the three species (inter and intraspecific cannibalism), while amphipods displayed predation of copepods, epilithic biofilm grazing, root shredding and sediment filter feeding.

Bulk stable isotope analysis

δ¹³C and δ¹⁵N SIA on bulk homogenised samples of sediment, roots and stygofauna (respectively 1.28 mg, 0.08–0.14 mg and 0.63–2.79 mg per samples, S2 Table) were performed at the Australian Nuclear Science and Technology Organisation (ANSTO, Sydney). Samples were loaded into tin capsules and analysed with a continuous flow isotope ratio mass spectrometer (CF-IRMS, Delta V Plus, Thermo Scientific Corporation, U.S.A.), interfaced with an elemental analyser (Thermo Fisher Flash 2000 HT EA, Thermo Electron Corporation, U.S.A.) following the procedure of Mazumder et al. [41]. δ¹³C values are reported in per mil (‰) relative to the Vienna Peedee Belemnite (VPDB), while δ¹⁵N values are reported relative to reference N₂ of known nitrogen isotopic composition (in ‰), previously calibrated against the AIR international isotope standard. δ¹³C POC (0.6 mg, S2 Table) was analysed at the Western Australian Biogeochemistry Centre at The University of Western Australia using a GasBench II coupled with a Delta XL Mass Spectrometer (Thermo-Fisher Scientific). Results have a precision of ± 0.10 ‰, and are reported relative to the NBS19 and NSB18 international carbonate standard [42].

Single amino acid carbon and nitrogen isotope analysis

δ¹³C CSIA

Stygofaunal samples (0.16–2.89 mg per sample, S2 Table) were hydrolysed under vacuum with 0.5 to 1 ml of amino acid-free 6 M HCl (Sigma-Aldrich) at 110 ˚C for 24 h. The protein hydrolysates were dried overnight in a rotary vacuum concentrator and stored in a freezer. Prior to analysis, the samples were resolved in Milli-Q water and 10 μl of 1-mmol solution of 2-aminoisobutyric acid (Sigma-Aldrich) was added as internal standard. The sample stock had a concentration of approximately 8 to 10 mg/ml, which was further diluted as needed. Single amino acid carbon isotope analysis was carried out at the La Trobe Institute for Molecular Sciences (LIMS, La Trobe University, Melbourne, Australia) using an Accela 600 pump connected to a Delta V Plus Isotope Ratio Mass Spectrometer via a Thermo Scientific LC Isolink (Thermo Scientific).

The amino acids were separated using a mixed mode (reverse phase/ion exchange) Primesep A column (2.1 x 250 mm, 100°C, 5 μm, SIELC Technologies) following the chromatographic method described in Mora et al. [43], after Smith et al. [44]. Mobile phases are those described in Mora et al. [45]. Percentage of Phases B and C in the conditioning run, as well as flow rate of the analytical run and timing of onset of 100% Phase C were adjusted as needed. Samples were injected onto the column in the 15 μl—partial loop or no waste—injection mode, and measured in duplicate or triplicate.

δ¹⁵N CSIA

CSIA nitrogen analyses were undertaken at the Organic Geochemistry Unit of the University of Bristol, UK. To extract the AAs, crushed samples (2.47–5.19 mg per sample, S2 Table) were hydrolysed in culture tubes (6 M HCl, 2 ml, 100°C, 24 h). A known quantity of norleucine (1 mg mL–1 in 0.1 M HCl) was added to each sample as an internal standard prior to hydrolysis. After heating, the tubes were allowed to cool then after centrifugation (3000 rpm, 5 min) the supernatant containing the hydrolysate from each tube was transferred to a clean culture tube and dried under N₂ whilst being heated to 70°C. Once dry, each sample was re-dissolved in 0.1 M HCl and stored in the dark at -18°C until required for analysis.

The derivatisation procedure followed Styring et al. [46] and included isopropylation, with a 4:1 mixture of 2-propanol and acetyl chloride heating to 100°C for 1 hour, the reaction was quenched by rapidly cooling in a freezer. After removing the residual solvents under N₂, acetylation of the amino group was achieved by adding a 5:2:1 mixture of acetone, triethylamine and acetic anhydride then heating to 60°C for 10 minutes before being allowed to cool. The derivatised AAs were isolated via liquid-liquid separation, residual solvent being removed by evaporating under N₂. Samples were again stored at -18°C until required for analysis.

A Thermo Finnigan Delta Plus XP isotope ratio mass spectrometer (Thermo Scientific, Bremen, Germany) was used to determine the δ¹⁵N values of derivatised AAs. The mass spectrometer (EI, 100 eV, three Faraday cup collectors for m/z 28, 29 and 30) was interfaced to a Trace 2000 gas chromatograph via a Combustion III interface (CuO/NiO/Pt oxidation reactor maintained at 980°C and reduction reactor of Cu wire maintained at 650°C), both from Thermo Scientific.

Samples were dissolved in ethyl acetate and 1μl of solution was injected via a PTV injector. Helium at a flow of 1.4 ml min^–1 was used as the carrier gas and the mass spectrometer source pressure was maintained at 9 X 10–4 Pa. The separation of the AAs was accomplished using a DB-35 capillary column (30 m X 0.32 mm i.d., 0.5 mm film thickness; Agilent Technologies, Winnersh, UK). The oven temperature of the GC started at 40°C where it was held for 5 min before heating at 15°C min^–1 to 120°C, at 3°C min^–1 to 180°C, at 1.5°C min^–1 to 210°C and finally at 5°C min^–1 to 270°C and held for 1 min. A Nafion dryer removed water and a cryogenic trap was employed to remove CO₂ from the oxidised and reduced analyte.

All the δ¹⁵N values are reported relative to reference N₂ of known nitrogen isotopic composition, previously calibrated against the AIR international isotope standard, introduced directly into the ion source in four pulses at the beginning and end of each run. Each reported value is a mean of duplicate δ¹⁵N determinations. A standard mixture of derivatised AAs of known δ¹⁵N values was analysed every three runs in order to monitor instrument performance.

Data treatment and statistical analysis

Only AAs that returned results for each taxon were considered. EAA and NEAA were separated according to the classification provided by Boudko [47]. EAAs were used in the interpretation of carbon flows—and potential shifts in OM incorporations—because they persist through the trophic chain [48] due to the little fractionation they undergo when incorporated into consumer’s tissue [49]. NEAA, which are subjected to much greater fractionation because of their de novo biosynthesis mainly from intermediates of the Krebs cycle (serine (Ser), glycine (Gly) and alanine (Ala)) and glycolysis (glutamic acid (Glx), aspartic acid (Asx) and proline (Pro) [50], were compared to EAA to investigate taxa-specific carbon isotopic trends (biosynthesis vs assimilation through diet) across the two rainfall periods (LR and HR).

All the statistical analyses were performed in R software version 3.6.0 (Development-Core-Team, 2016). Analysis of variance (ANOVA) coupled with Tukey’s HSD pairwise comparisons (R-package ‘stats’) were employed to inspect significant differences between δ¹³C_EAA (Val, Phe and Arg) and δ¹³C_NEAA (Krebs (Ser, Gly and Ala) and glycolysis (Asx, Glx and Pro) cycles) within the different rainfall conditions (LR and HR). Principal component analyses (PCA, R-package ‘vegan’) and Linear Discriminant Analysis (LDA, R-package ‘vegan’) among EAA was performed to explore sample distribution in the multi-dimensional space. Determination of EAA driving sample variability in the PCA was carried out via function fviz_contrib (R-package ‘factoextra’).

Trophic positions (TP) were calculated using the methodology reported by Chikaraishi et al. [33]:

Incorporation of source carbon from terrestrial vegetation has previously been reported at Sturt Meadows, with roots from surficial saltbush vegetation (C3 metabolism) frequently found in the groundwater [40]. β was accordingly assigned the value of +8.4 ± 1.6 ‰, which is the established value for aquatic food webs involving C3 plants [31]. Although other carbon sources are possible in groundwaters, as they are not established in this system, a conservative approach has been taken in using the value of an evidenced source. TDF was assigned the value of 7.6 ± 1.2‰, based on Steffan et al. [51] who showed it did not vary across trophic levels one to four in multiple controlled-feeding experiments, and for trophic levels one to five in a natural food chain, using terrestrial arthropod species [28].

Pairwise comparisons for δ¹⁵N were carried out with the same approach as for the carbon CSIA data. Robustness and consistency between CSIA and SIA data from beetles and amphipods were inspected using Pearson correlations (function rcorr in R-package ‘Hmisc’). SIMM (Stable Isotope Mixing Models, R-package ‘simmr’) were then applied to establish dietary proportions of the key ecological taxa (Fig 2). Since a specific trophic discrimination factor has not been calculated for stygofauna, we used the widely accepted values of 3.4 ± 2 ‰ for nitrogen and 0.5 ± 1 ‰ for carbon [52]. Markov chain Monte Carlo (MCMC) algorithms were used for simulating posterior distributions in SIMM, and MCMC convergence was evaluated using the Gelman-Rubin diagnostic by using 1.1 as a threshold value for analysis validation.

Results

Stygofaunal carbon fluxes

During LR, δ¹³C average values of AAs (δ¹³C_NEAA[LR] and δ¹³C_EAA[LR]) spanned from -31.52 ‰ (Phe) to -5.72 ‰ (Gly). Similar values were found under HR conditions (δ¹³C_NEAA[HR] and δ¹³C_EAA[HR]), ranging from -31.55 ‰ (Phe) to -4.92 ‰ (Ser) (Table 1).

δ¹³C_NEAA[LR] average values varied from -26.36 ‰ (Ala) to -5.72 ‰ (Gly), similar values to δ¹³C_NEAA[HR] spanning from -26.28 ‰ (Ala) to -4.92 ‰ (Ser). Overall, δ¹³C_EAA showed trends towards more negative values than δ¹³C_NEAA, which is involved in both Krebs and glycolytic cycles (ANOVA, P <0.005). This is consistent with the enrichment of NEAA during biosynthesis in the organism. With the exception of water mites (OR), P. macrosturtensis larvae (Blv), Y. sturtensis (AM2), S. bradfordae (AM3) under LR, and oligochaetes (OL) under HR, pairwise comparisons between δ¹³C_EAA and δ¹³C_NEAA confirmed a shift towards more negative values across the stygofaunal community (Table 1).

Neither PCAs nor LDAs on EAA distinguished different clusters within taxa or main groups (adult and larval beetles, and amphipods) nor among different rainfall periods (LR and HR). All three EAA correlated positively and significantly (P < 0.005), with phenylalanine and valine being the most informative AAs explaining the isotopic variability across stygofauna (~70%). δ¹³C values of valine (δ¹³C_Val) and phenylalanine (δ¹³C_Phe) show that, with the exception of P. microsturtensis (S) and S. bradfordae (AM3), the entire stygofaunal community experienced a significant shift towards more ¹³C-depleted values under HR (MANOVA, P < 0.005) (Fig 3, Table 2).

Within the significant trends, P. macrosturtensis adults and larvae (B and Blv) showed the smallest change in carbon values (B: δ¹³C_Val+Phe = -2.91; Blv: δ¹³C_Val+Phe = -3.89) between rainfall regimes, while amphipods S. axfordi and Y. sturtensis (AM1 and AM2) showed the largest depletion (AM1: δ¹³C_Val+Phe = -12.59 ‰; AM2: δ¹³C_Val+Phe = -10.10 ‰), suggesting differential carbon incorporations under HR conditions.

δ¹⁵N and trophic levels

δ¹⁵N_Glu average values varied between 15.4±0.4‰ (AM3_[HR]) and 22.31±0.29‰ (M_[HR]), while δ¹⁵N_Phe values ranged from 10.67±0.45‰ (AM3_[HR]) to 14.53±0.06‰ (M_[HR]). When converted to trophic positions, the stygofaunal community at Sturt Meadows shows a truncated trophic chain, clustering around the secondary consumer level (Fig 4).

Under LR conditions, P. macrosturtensis larvae (Blv) show the highest trophic position (TP = 3.33±0.02), while water mites (OR) sit at the lowest (2.78±0.09). Under HR conditions, P. microsturtensis adults (S) have the highest trophic position (3.27±0.01), whilst S. bradfordae (AM3) show the lowest value (TP = 2.73±0.01). Due to the low abundances it wasn’t possible to analyse biochemical fingerprints from water mites (OR_[HR]: 37 individuals) and P. microsturtensis larvae (Slv_[HR]: 10 individuals) during the wet season (HR) (S1 Table).

Overall, adult beetles (B, M and S) revealed higher trophic levels (TP>3) than amphipods (AM1, AM2 and AM3, TP<3). However, B, M and S did not show statistically higher values than AM1 under LR, the same pattern seen in P. mesosturtensis (M) under HR.

S. bradfordae (AM3) and P. macrosturtensis larvae (Blv) are the only organisms to show a statistically significant change in their TP values between LR and HR (Table 3), with both decreasing trends.

Food web dynamics

CSIA-based TP correlated significantly with SIA δ¹⁵N and δ¹³C values both under LR (P<0.05 in both cases) and HR conditions (P<0.01 and P<0.05 respectively). Under the latter conditions, δ¹³C_Val values correlated significantly with CSIA-based TP (P<0.05). Copepods are generally thought to sit at the base of the food web [53,54]. However, these were analysed only via bulk SIA due to organism and sample size, and so could not be included in the TP analysis. They showed more ¹³C-depleted δ¹³C (cyclopoids: δ¹³C_LR = -20.5‰, δ¹³C_HR = -21.9‰; harpacticoids: δ¹³C_LR = -20.6‰, δ¹³C_HR = -23.5‰) and enriched δ¹⁵N (cyclopoids: δ¹⁵N_LR = 13.9‰, δ¹⁵N_HR = 14.5‰; harpacticoids: δ¹⁵N_LR = 11.9‰, δ¹⁵N_HR = 15.8‰) values under HR (Fig 5A and 5B), indicating that the change in rainfall regimes could play a role in stygobiotic meiofaunal biochemical incorporations.

Amphipods AM1 and AM2 sat at the base of the trophic web under both rainfall conditions (TPs always below three, Table 3), and SIA carbon values (δ¹³C) confirmed a shift–already pinpointed via CSIA—towards more ¹³C-depleted carbon sources under HR. AM3, the smallest and rarest amphipod species in the calcrete, did not allow bulk SIA analyses due to the low abundances detected (LR (average value between LR1 and LR2): 27; HR: 19, S1 Table).

With respect to dietary preferences, for the amphipod S. axfordi (AM1), mixing models suggest that roots (and hosted microbial flora) contributed the greatest proportion (50%) during low rainfall conditions (LR) (Fig 6). The remaining diet was composed of POC (16.1%, derived from allochthonous carbon incorporations, a potential organic source for microbes), copepods (harpacticoids (13.9%) and cyclopoids (11.9%)) and sediment (8.1%) (i.e. OM laying at the bottom of the aquifer or epilithic biofilms).

Under HR conditions, the POC dietary contribution reached 66.1%, while roots plummeted to 3.3% (Fig 6). Overall, amphipod Y. sturtensis (AM2) showed the same dietary patterns as AM1.

Adult beetles P. macrosturtensis (B) and P. mesosturtensis (M) show only slight depletions in their isotopic values during HR in bulk δ¹³C and δ¹⁵N SIA, in contrast to the larger changes seen in the CSIA data. P. microsturtensis (S), which showed an isotopic enrichment in CSIA, counter to the rest of the community, shows similar behaviour to P. macrosturtensis (B) and P. mesosturtensis (M) in the SIA (S3 Table). All the three species show similar dietary preferences in mixing models across the rainfall periods (S4 Table). While diets were dominated by amphipods AM1 and AM2 during the LR period (B: 39.9%, M: 49.3% and S: 47.9% (Fig 7)), predation/scavenging of sister beetle species accounted for the biggest dietary proportions during the wet season (B: 52.9%; M: 49.4%; S: 41.9% (Fig 7)).

Mixing models indicate that P. macrosturtensis larvae (Blv), which showed the biggest shift in trophic position, has a preference for amphipods S. axfordi (AM1) and Y. sturtensis (AM2) under LR conditions (accounting for 52% of the diet contributions), but also consumes a range of other organisms (Fig 8). During HR, Blv’s diet is dominated by the two amphipod species, accounting for 79.6% of food sources (Fig 8). Overall, these results indicate changes in amphipods (AM1 and AM2) diet preferences linked with different OM inputs, coupled with enhanced species-specific predatory pressures from Blv under HR conditions.

Discussion

Shifts in basal OM assimilation

δ¹³C_EAA data suggest that the stygofaunal community at Sturt Meadows experienced a seasonal shift in carbon flows during the wet season (HR). The overall tendency towards more ¹³C-depleted δ¹³C_Val and δ¹³C_Phe values indicate general stygofaunal discrimination against ¹³C sources in that season.

Three groups (OL, AM3 and S) showed a counter trend of ¹³C enrichment in the EAAs during HR. Of these, the easiest to account for are the oligochaetes (OL) which also showed increased abundances (χ² = 6.7698, P < 0.05) of individuals ranging from 2 (LR1) and 1 (LR2) to 21(HR) (S1 Table), indicating ideal conditions for the taxon during the wet season. As detritovores, oligocheates may be expected to preferentially consume more degraded, and so ¹³C-enriched, OM. The enrichment in S. bradfordae (AM3) and P. microsturtensis (S) is harder to explain at this stage. The low abundance of S. bradfordae (AM3) means that, like the oligochaetes, it was not included in the SIA analysis and less data is available. This taxon, together with P. macrosturtensis larvae (Blv), Y. sturtensis (AM2) and water mites (OR), lacked statistically significant differences when δ¹³C values of EAA during LR are compared with those of NEAA involved in the glycolytic cycle. Newsome et al. [50] indicated that δ¹³C values of NEAA from diets of omnivorous animals reflect de novo synthesis but also dietary incorporations. Differential routing of macromolecules by consumers [55] are one possible contributor to our results. However, to date isotopic routing hypotheses have been tested only in vertebrates [30,56], with the study of metabolic pathways in aquatic invertebrates largely unexplored. Further CSIA investigations involving species-specific bio-assimilation processes within the stygofaunal community are needed to provide a more accurate understanding of the biochemical dynamics regulating this system.

In line with our general data trends, Hartland et al. [15], who reported consistent depletion in δ¹³C stygofaunal values within OM-enriched groundwaters via sewage contamination, concluded that stress-subsidy gradients in groundwaters trigger profound changes in stygofaunal assemblages and have the potential to trigger shifts in feeding habits. Rainfall events trigger OM inflows which constitute high quality carbon sources for aquatic biota in groundwaters [16, 57].

Reiss et al. [26] demonstrated a strong link between nutrient inputs (mainly DOC) and groundwater microbe functional and metabolic richness after a major flooding event. Unfortunately, their methodology did not allow for corresponding macrofaunal trends to be identified. Nonetheless, microbially-derived OM incorporation by stygofauna has been reported in a number of groundwater ecology studies [16,21,58], and the biochemical importance of this linkage is widely accepted.

A key role in the observed trends at Sturt Meadows is played by amphipods which, together with copepods, are recognised as crucial actors in transferring OM to the upper stygofaunal trophic levels [16]. Specialized trophic habits in amphipods include epigean predation [10], detritivory [59], parasitism [60], biofilm grazing [61] and necrophagy [62]. Several studies have reported high degrees of trophic opportunism [54] and plasticity [63], allowing amphipods crucial shifts in feeding modes. Concurrently, niche partitioning has been addressed as a key mechanism to reduce intraspecific competition in ecosystems shaped by scarce nutrient availability [64]. However, our results do not show any conclusive evidence of epigean amphipod niche partitioning, with amphipods S. axfordi (AM1) and Y. sturtensis (AM2) showing the same dietary patterns. Overall, the isotopic data support the concept of opportunistic behaviours linked with changes in resource availability as a result of different rainfall regimes.

The HR event triggered substantial changes in the dietary proportions of S. axfordi (AM1) and Y. sturtensis (AM2), with notable decreases in root input and increases in POC. The extent of direct plant matter consumption by stygobionts reported in the literature–particularly by amphipods, which are facultative shredders–is both site and species-specific. Jasinska et al. [65] found that aquatic root mats were a key food source for a biodiverse cave fauna hosted by a shallow groundwater stream in Western Australia. Conversely, Navel et al. [66], reported the widely distributed amphipod species Niphargus rhenorhodanensis having preferential OM collector/gatherer feeding habits. In another study, Simon et al. [67] suggested that wood inputs played a role as indirect source of OM consumed by the ephilitic microbial mats which were ultimately targeted by common Gammarus amphipods.

At Sturt Meadows, a plausible explanation for the patterns observed is that during the dry season epigean amphipods rely on a more omnivorous diets where roots falling from the surface, and associated microbial and fungal biota, provide a substantial food source. Conversely, the wet season triggers inflows of replenished carbon (¹³C-depleted POC) that might fuel biological turnovers in microbiological activity, and POC-attached microflora may be ultimately targeted by epigean amphipods. These assumptions are in line with the finding reported by Brankvotis et al. [16], and support the concept that grazers play a crucial role in sustaining the functional diversity in groundwaters. The importance of plant matter input during at least part of the year is supported by a previous bulk SIA investigation at Sturt Meadows [40], which also suggested that terrestrial sources of carbon, mainly DOC, reached the aquifer via percolation and play a crucial role in energy flows within the system. It is worth noting that our δ¹³C values of decarbonated sedimentary fractions (referred above as ‘sediment’) were less ¹³C-depleted than those in other groundwater investigations ([68,69,70], and had ranges close those for dissolved organic carbon (DIC) in the region ([71,72,73]. Portillo et al. [74] reported karst microbial growth induced by both carbonate precipitation and dissolution, suggesting the inclusion of inorganic carbon within the estimation of global carbon budgets in groundwaters. In line with this work, Chapelle [75] reported in situ DIC production as a result of microbial metabolism involved in the dissolution of carbonate material in the black Creek aquifer (California, USA). Our results suggest that carbonate assimilation and/or dissolution processes are likely to occur in the sedimentary deposits of the aquifer, transferring an inorganic carbon isotopic fingerprint into the decarbonated and organic fractions. This can be tested in future by further functional studies on the microbial community [76] at the site.

Copepods (C and H) showed high δ¹⁵N values compared to the rest of the stygofaunal community (S3 Table), suggesting alternative nitrogen sources linked to different microbial baselines. Copepods act as energy drivers in recycling nitrogen via ingestion of sediment and attached bacteria [12], with ammonia (NH₃), together with nitrates (NO₃^-), being an essential nutrient and energy source for subterranean microorganisms [77]. At Sturt Meadows aquifer, where ammonia levels are considerably higher than the natural concentrations [38], proliferation of selectively grazed ammonia-oxidising bacteria (AOB) might have played a key role in triggering the enriched δ¹⁵N values in copepods.

The present study is constrained by its focus on stygofauna and therefore cannot provide direct evidence of microbially-derived ecological shifts. Future research needs to combine stygofaunal and microbial investigations to create a complete picture of the ecosystem. CSIA and functional genetic studies on microbes and copepods would also help define transitions from microflora to stygofauna, a process that has been understudied so far. Recent promising investigations in surface terrestrial [34] and aquatic [78] environments suggest a design—carbon fingerprinting—based on the incorporation of isotopic data into multifactorial mixing models that allow specific elucidation of bacterial sources in diets. Overall, despite the methodological challenges posed by groundwaters, isotopic data on stygofaunal carbon fluxes provides baseline knowledge that help untangle the intricate biochemical dynamics regulating subsurface environments.

Trophic interactions

Our data on nitrogen CSIA pinpointed two main trophic levels marked by a small but clear separation between the top predators—adult beetles (B, M and S)—and primary consumer amphipods (AM1, AM2 and AM3) under both rainfall conditions. Compared to other ecosystems [29], the Sturt Meadows aquifer shows a very simple and truncated trophic web dominated by omnivorous habits. This is consistent with previous assumptions [67] due to the lack of primary producers [9] and scarce nutrient availability [79].

Within subterranean beetles, the smallest species P. microsturtensis (S), together with P. macrosturtensis (B), sit at the top of the trophic chain during HR (Table 3). Under those conditions, increased oxygen levels [38] may play a role in shaping changes in stygofaunal niche occupation. Subterranean beetles’ body size has been found to drive differential physiological responses to increased exoskeleton respiration rates (inversely proportional to the body size) which ultimately affect the ability to allocate energy for breeding and foraging [80]. As the smallest species P. microsturtensis (S) can adapt their metabolism more quickly than direct competitor sister species P. mesosturtensis (M) under favourable conditions—such as HR regimes–they are more likely to show shifts in ecological niche occupation [38]. This trend, combined with the group feeding tendency of P. microsturtensis (S) beetles [40], indicates higher efficiency in activating more intensified predatory strategies when compared to P. mesosturtensis (M).

Dytiscidae beetle larval stages—commonly referred as ‘water tigers’—are ferocious carnivores [81] with extremely opportunistic feeding behaviours involving scavenging and cannibalism [82]. At Sturt meadows, the third instar of blind P. macrosturtensis larvae (Blv) has a considerably bigger head capsule—paired with elongated mouthparts—than adult stages (Fig 9). These morphological features are likely to provide ethological advantages for non-visual predacious habits within light-less environments such as groundwaters [83]. This is consistent with stable isotope data from LR conditions positioning Blv at the top of the trophic web (Table 3).

Overall, modelled dietary contributions of P. macrosturtensis larvae (Blv) indicated a preference for amphipods (AM1 and AM2) coupled with residual cannibalism/scavenging (B, M and S) and predation of copepods (C and H) (Fig 8). Under HR, P. macrosturtensis larvae (Blv) showed the biggest drop in TP compared to LR (TP_LR-HR = -0.3), which can be explained by an increased predatory focus on amphipods, and reduced feeding on secondary consumer sister species.

Previous work on surficial Dytiscidae larval stages published by Inoda et al. [84] stressed the importance of prey recognition through smell. According to their results, prey density was found not to be a key factor in shaping feeding behaviours, and self-other recognition played a role in group feeding. Overall, these findings indicated prevention of cannibalism through scent recognition. In groundwater, with total darkness and high influence of OM inputs on population dynamics [85,86], these patterns are likely to be strengthened. We suggest that the shifts in Blv predation seen in our results are dictated by a combination of chemical recognition and increased likelihood of encountering prey (amphipods) driven by enhanced resource availability (OM) during HR periods.

The role of bottom-up vs top-down forces in natural communities has been a cornerstone issue in the field of trophic ecology since the first empirical investigations [87]. Despite the controversy generated by the debate, there is now consensus that both forces act simultaneously on populations. This reinforces the need for whole system studies considering the interaction between heterogeneous (biotic and abiotic) forces and their effect on communities [88,89,90].

Our study, in line with a number of other investigations in the field [26,91] confirms that rainfall events via water advection are key drivers in defining energy flows and ecological patterns within resource-limited environments. We suggest that OM-driven bottom-up regulations, increasingly accepted as driving factors shaping population dynamics in aquifers [92], shaped the shifts in feeding behaviours among amphipod taxa in the calcrete. However, despite the beneficial conditions for primary consumers triggered by increased nutrient availability (i.e. microbial biofilms) and better environmental settings (i.e. increased oxygen, [38]), a decrease in amphipod populations under HR indicates the existence of additional ecological factors.

Top-down forces (i.e. natural predators), widely studied in surface aquatic ecosystems [93,94], have hardly been addressed in groundwater. Previous genetic investigations at Sturt Meadows pinpointed predatory pressures from beetles on amphipods and copepods [95] and reported a lack of trophic niche partitioning among the Paroster species. In another study, Hyde [96] reported evidence from metagenomics data suggesting that subterranean blind beetles at Sturt Meadows feed on both prey invertebrates and their sister species. Our isotope results support these hypotheses, indicating opportunistic predaceous habits in the calcrete, mixed with scavenging/cannibalism. However, substantial uncertainty remains about the magnitude of interspecific predatory pressures among Paroster sister species, and further species-specific lab experiments are needed to investigate these ethological aspects.

Biochemical functional role interpretation coupled with abundance data suggests that bottom-up population dynamics are counterbalanced in the system by top-down forces. Increased numbers of top predators (adult beetles B, M and especially S) were paired with a decrease of key prey items (amphipods AM1, AM2 and AM3) when HR is compared with the dry season (LR1 and LR2) (Fig 10).

In light of the dynamics shown by our isotope data, we suggest that the reported shift in amphipods (AM1 and AM2) carbon incorporation during HR might have triggered changes in their ecological behaviour, exposing them to increased predatory pressures from the top predator Paroster beetles (B, M and S). However, given the high degree of opportunistic behaviour reported by stygofauna [14], further investigations on species-specific ethological dynamics would be helpful to infer community dynamics.

The number of third instar dytiscidae larvae ‘Blv’ did not vary across sampling campaigns, suggesting differential ecological niche occupations. Previous investigations on Paroster larvae detected three instars before pupating, with the first two occupying a reduced proportion of their lifetime [83]. Future investigation of early stages of larval developments are needed to establish if potential population blooms (i.e. mass reproduction) are linked with contrasting recharge periods.

Conclusions

The application of CSIA and SIA allows elucidation of the trophic dynamics shaping stygofaunal communities in an arid zone calcrete aquifer. Rainfall acts as a driver in regulating both top down and bottom up changes in dietary habit. Subterranean invertebrate population dynamics are notoriously hard to investigate due to sampling obstacles and a current lack of knowledge around stygofaunal biological cycles [7,36]. However, our isotopic results allow a greater insight into the food web dynamics and the biogeochemical forces that shape them than has previously been possible. Further investigations involving higher numbers of samples from more biodiverse systems or complex trophic assemblages (i.e. alluvial aquifers) will help refine the approach. The incorporation of qualitative analyses such as DNA metabarcoding would also complement quantitative isotopic methods to provide crucial insights into processes (i.e. cannibalism) and key driving forces (i.e. niche partitioning) that are hard to detect via one method alone. Lastly, investigation of nitrogen sources and their isotopic changes would open up the nitrogen data collected to interpretation beyond trophic position.

Groundwater environments are fundamentally important to ecosystems, communities and industry, and a robust understanding of their ecosystem dynamics is essential to accurately assess environmental impacts, whether anthropogenic, or climatic. Isotopic data, especially if combined in multidisciplinary studies with other parameters [22] has a key role to play in elucidating previously hard to investigate function within these cryptic systems.

Supporting information

S1 Table [xlsx]
Abundance data of the different stygofaunal taxa at Sturt Meadows detected during the sampling campaigns LR1, LR2 and HR.

S2 Table [xlsx]
Sample weights (mg) for SIA and CSIA analyses.

S3 Table [xlsx]
δ13C and δ15N values of stygofauna, roots, sediment and POC during LR and HR.

S4 Table [xlsx]
Dietary proportions of (B), (M) and (S) under LR and HR conditions.