Rainfall as a trigger of ecological cascade effects in an Australian groundwater ecosystem

Study area

The field work was carried out at the Sturt Meadows 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, see Fig. 1a). The study area is a calcrete aquifer lying in the Raeside paleodrainages in the Yilgarn region of Western Australia (Fig. 1a). The vegetation of the area is Acacia woodlands, primarily Acacia aneura (F.Muell. ex Benth.), and is subjected to combined grazing pressure from domestic stock, feral animals and macropods. The aquifer is accessible through a bore grid comprising 115 bore holes of between 5 and 11 m depth (Fig. 1b).

Three sampling campaigns were carried out, two of them (LR1: 26/07/2017 and LR2: 07/11/2017) corresponding to low rainfall periods²⁴ and one during the wet season (high rainfall, HR; two consecutive days of sampling collection: 17-18/03/2018) (Supplementary Fig. 1). The well-studied stygofaunal community of the area is composed of 11 main stygofaunal taxa belonging to five Classes: Oligochaeta (family Tubificidae (Vejdovský 1884)), subcohort Hydrachnidia, Maxillopoda (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)), Malacostraca: Amphipoda (species Scutachiltonia axfordi (King, 2012), Yilgarniella sturtensis (King, 2012) and Stygochiltonia bradfordae (King, 2012)) and Insecta: Coleoptera: Dytiscidae (species Paroster macrosturtensis (Watts & Humphreys 2006), Paroster mesosturtensis (Watts & Humphreys 2006) and Paroster microsturtensis (Watts & Humphreys 2006) and respective larvae).

Field work procedures and sample preparation

Given the sensitivity of the hydrological dynamics in shallow calretes^23,25, extensive water extraction along the bores was avoided, and preliminary tests on the bores with the highest water depth were carried out to quantify potential risk of dewatering the calcrete. During the field campaigns LR2 and HR, 20 water samples in total (two samples for stable isotope analysis on DOC (Dissolved Organic Carbon) and DIC (Dissolved Inorganic Carbon), three samples for radiocarbon analysis on DOC, one sample for radiocarbon analysis on DIC, and two samples for stable isotope and radiocarbon analyses on POC (Particulate Organic Carbon)) were collected from bores D13 and W4 (Fig. 1b), which are representative of the two main geological conformations of the area—calcrete (W4) and clay (D13) (Supplementary Fig. 2)—and host stable hydrological and biotic conditions⁷. Water samples were collected using a submersible centrifugal pump (GEOSub 12 V Purging Pump) after wells were purged of three well-volumes and stabilisation of in-field parameters was observed, according to the methodology in Bryan et al.²⁶.

Samples for ¹⁴C_DIC analysis were filtered through 0.45 μm filters and collected in 1 L high density poly-ethylene (HDPE) bottles. δ¹³C_DIC samples were filtered through 0.2 μm filters, collected in 12 mL glass vials (Exetainers) and refrigerated after sampling. δ ¹³C_DOC samples were filtered through 0.2 μm filters, collected in 60 mL HDPE bottles and frozen after sampling.¹⁴C_DOC samples were filtered through 0.2 μm filters, collected in 3 L HDPE bottles and frozen after sampling.

In order to investigate ¹⁴C and δ ¹³C content of POC, two extra liters were collected from the same bores (D13, W4) and kept frozen (− 20 °C) until further analyses. ¹⁴C_POC δ¹³C_POC samples were then filtered through pre-combusted GF/F filters (12 h at 450 °C), washed with 1.2 N HCl to remove any inorganic carbon, and subsequently dried at 60 °C for 24 h. All samples were closed with sealing tape after collection to limit atmospheric exchange and kept in darkness.

Temperature, pH, ORP, salinity, DO and depth were measured in situ using a portable Hydrolab Quanta Multi‐Probe Meter across 30 bores during LR1, LR2 and HR²³ (presented in Supplementary Table 8). Adult and larval stygofaunal specimens were collected from the same bores by hauling a weighted plankton net (mesh 100 μm²⁷) five times through the water column (Fig. 1b). All biological samples were kept frozen (− 20 °C) in darkness until laboratory processing. Individual organisms were counted and identified (and consequently separated) to the lowest taxonomic level via optical microscopy and reference to specific taxonomic keys. Plant material, sediment samples and fauna were each separated during the sorting in the laboratory and each taxon pooled according to sampling campaign (LR1, LR2 or HR) and subsequently washed with Milli-Q water to remove surface impurities from their bodies. Sediment samples were soaked in acid (0.1 N HCl) to remove inorganic carbon, and together with the other samples were then oven dried at 60 °C overnight and ground until obtaining a homogeneous fine powder and stored at − 20 °C until further analyses.

Previous investigations on the ecological niche trends at Sturt Meadows indicated that all stygofauna characterize similar niche occupations under low rainfall regimes (LR1 and LR2)²³. Stygofaunal specimens from the two low rainfall sampling events were combined to form sample LR to address the competing requirements between isotopic detection limits, analytical replicates and cost, while maintaining the main taxonomic and biological classifications⁷.

Bulk isotope and ¹⁴C analyses

Water δ¹³C_DIC and δ¹³C_POC isotopic ratios were analysed by Isotope Ratio Mass Spectrometer—WABC at The University of Western Australia using a GasBench II coupled with a Delta XL Mass Spectrometer (Thermo-Fisher Scientific)—and the results, with a precision of ± 0.10‰, were reported as ‰ deviation from the NBS19 and NSB18 international carbonate standard²⁸.

δ¹³C_DOC isotopic ratios of waters were analysed via Liquid Chromatography Isotope Ratio Mass Spectrometer (LC-IRMS) at the La Trobe Institute for Molecular Sciences (LIMS, La Trobe University, Melbourne, Australia) comprising an Accela 600 pump connected to a Delta V Plus Isotope Ratio Mass Spectrometer via a Thermo Scientific LC Isolink (Thermo Scientific).

C and N bulk SIA on homogenised samples of sediment, roots, stygofauna and copepods (cyclopoids and harpacticoids) were performed at the Australian Nuclear Science and Technology Organisation (ANSTO, Sydney, Australia). Samples were loaded into tin capsules and analysed with a continuous flow isotope ratio mass spectrometer (CF-IRMS), model Delta V Plus (Thermo Scientific Corporation, U.S.A.), interfaced with an elemental analyser (Thermo Fisher Flash 2000 HT EA, Thermo Electron Corporation, USA) following the procedure published by Mazumder et al.²⁹.

For radiocarbon analyses, samples (sediment, roots, copepods, ants, stygofauna) were treated with 1 M HCl for 2 h to remove all possible carbonate contamination. These pre-treated samples together with ¹⁴C_POC, ¹⁴C_DOC and ¹⁴C_DIC samples were subjected to CO₂ extraction and graphitization following the methodology published by Hua et al.³⁰. ¹⁴C content of samples was determined by means of the Accelerator Mass Spectrometry (AMS) at ANSTO.

Carbon CSIA

Carbon CSIA followed the procedure described in Saccò et al.⁷. Samples of roots and stygofaunal specimens 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 dissolved in Milli-Q water and 10 μL of 1-mmol solution of 2-aminoisobutyric acid (Sigma-Aldrich) 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 × 250 mm, 100 °C, 5 μm, SIELC Technologies) following the chromatographic method described in Mora et al.³¹ after Smith et al.³². Mobile phases are those described in Mora et al.³³. 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.

To elucidate carbon flows through the stygofaunal community we focused on the essential amino acids Valine (Val), Phenylalanine (Phe) and Arginine (Arg), as these compounds must be integrated through diet and cannot be synthetised internally by the fauna^14,34. In addition, to distinguish between terrestrial and aquatic carbon sources, the ratio between Val and Phe signals (δ¹³C_Val-Phe), a widely employed index in archaeology and freshwater biology³⁵, was calculated for roots, water mites, aquatic worms, amphipods and beetles (larvae and adults).

Microbial taxonomic and functional gene analyses

Consumer amphipods (Scutachiltonia axfordi (AM1), Yilgarniella sturtensis (AM2), S. bradfordae (AM3)), cyclopoids and harpacticoids, together with predator stygobiotic beetles (Paroster macrosturtensis (B), P. mesosturtensis (M) and P. microsturtensis (S)) (see Saccò el at ⁷. for further details on the trophic characterisation of the stygofaunal community at Sturt Meadows), were used for gut microbiome bacterial 16S metabarcoding and microbial functional analysis. A total of 16 AM1, 16 AM2, 16 AM3, 20 cyclopoids and 20 harpaticoids and 20 of each one of the three Paroster species (B, M and S), were sorted into duplicates of stygobiotic pools of 3–5 individuals from both LR and HR events for DNA extraction. Prior to DNA extraction stygobitic animals (3–5 individuals per pool; n = 40) were placed in a petri dish containing ultrapure water and UV sterilized in a UV oven for 10 min to eliminate any bacterial species that may be contained on the exoskeleton as this study targeted the gut microbiome. Immediately post-UV treatment, the animals were placed into tissue lysis tubes with 180 μL tissue lysis buffer (ATL) and 20 μL Proteinase K and homogenised using Minilys tissue homogeniser (ThermoFisher Scientific, Australia) at high speed for 30 s. Lysis tubes, inclusive of two laboratory controls, were incubated at 56 °C using an agitating heat block (Eppendorf ThermoStat C, VWR, Australia) for 5 h.

Following the incubation, the analytical procedure was adapted from Saccò et al.²² and DNA extraction was carried out using DNeasy Blood and Tissue Kit (Qiagen; Venlo, Netherlands) and eluted off the silica column in 30–50 μL AE buffer. The quality and quantity of DNA extracted from each stygobitic pool was measured using quantitative PCR (qPCR), targeting the bacterial 16S gene. PCR reactions were used to assess the quality and quantity of the DNA target of interest via qPCR (Applied Biosystems [ABI], USA) in 25 μL reaction volumes consisting of 2 mM MgCl₂ (Fisher Biotec, Australia), 1 × PCR Gold Buffer (Fisher Biotec, Australia), 0.4 μM dNTPs (Astral Scientific, Australia), 0.1 mg bovine serum albumin (Fisher Biotec, Australia), 0.4 μM of each primer (Bact16S_515F and Bact16S_806R^36,37), and 0.2 μL of AmpliTaq Gold (AmpliTaq Gold, ABI, USA), and 2 μL of template DNA (Neat, 1/10, 1/100 dilutions). The cycling conditions were: initial denaturation at 95 °C for 5 min, followed by 40 cycles of 95 °C for 30 s, 52 °C for 30 s, 72 °C for 30 s, and a final extension at 72 °C for 10 min.

DNA extracts that successfully yielded DNA of sufficient quality, free of inhibition, as determined by the initial qPCR screen (detailed above), were assigned a unique 6–8 bp multiplex identifier tag (MID-tag) with the bacterial 16S primer set. Independent MID-tag qPCR for each stygobiotic pool were carried out in 25 μL reactions containing 1 × PCR Gold Buffer, 2.5 mM MgCl₂, 0.4 mg/mL BSA, 0.25 mM of each dNTP, 0.4 μM of each primer, 0.2 μL AmpliTaq Gold and 2–4 μL of DNA as determined by the initial qPCR screen. The cycling conditions for qPCR using the MID-tag primer sets were as described above. MID-tag PCR amplicons were generated in duplicate and the library was pooled in equimolar ratio post-PCR for DNA sequencing. The final library was size selected (160–600 bp) using Pippin Prep (Sage Sciences, USA) to remove any MID-tag primer-dimer products that may have formed during amplification. The final library concentration was determined using a QuBitTM 4 Fluorometer (Thermofischer, Australia) and sequenced using a 300 cycle V2 kit on an Illumina MiSeq platform (Illumina, USA).

MID-tag bacterial 16S sequence reads obtained from the MiSeq were sorted (filtered) back to the stygobitic pool based on the MID-tags assigned to each DNA extract using Geneious v10.2.5³⁸. MID-tag and primer sequences were trimmed from the sequence reads allowing for no mismatch in length or base composition.

Filtered reads were then input into a containerised workflow comprising USEARCH³⁹ and BLASTN⁴⁰, which was run on a high-throughput HPC cluster at Pawsey supercomputing centre. The fastx-uniques, unoise3 (with minimum abundance of 8) and otutab commands of USEARCH were applied to generate unique sequences, ZOTUs (zero-radius Operational Taxonomic Units) and abundance table, respectively. The ZOTUs were compared against the nucleotide database using the following parameters in BLASTN: perc_identity ≥ 94, evalue ≤ 1e−3, best_hit_score_edge 0.05, best_hit_overhang 0.25, qcov_hsp_perc 100, max_target_seqs = 5. An in-house Python script was used to assign the ZOTUs to their lowest common ancestor (LCA)⁴¹. The threshold for dropping a taxonomic assignment to LCA was set to perc_identity ≥ 96 and the difference between the % of identity of the two hits when their query coverage is equal was set to 1. Results on the microbial taxonomic diversity detected in ground water samples from a previous study on carbon inputs in the aquifer²² were incorporated in this work to allow qualitative comparison with the stygofaunal gut diversity.

To investigate functional activity involved in carbon cycling, the 16S metabarcoding data were fed to the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2 (PICRUSt2) software package to generate predicted metagenome profiles⁴². These profiles were clustered into Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthologs (KOs)⁴³ and MetaCyc pathway abundances⁴⁴ focusing on the relative abundances of four carbon metabolisms: carbon fixation in prokaryotes, carbohydrates, lipids and amino acid metabolisms. Relative abundance of pathways linked with methane, nitrogen and sulfur metabolisms were also investigated.

Statistical analyses

The Phyloseq package in R^45,46 was used to plot the ZOTU abundance for the styfofaunal specimens at the order level under LR and HR periods. The Statistical Analysis of Metagenomic Profiles (STAMP) bioinformatics software package was used to carry out Principal Components Analysis (PCA) to assess the differences between functional genomic fingerprints based on the KEGG orthologs function prediction between copepods (C and H) and amphipods (AM1, AM2 and AM3), and determine statistically significant results from the PICRUSt2 output⁴⁷. Clustering patterns according to rainfall periods (LR and HR) and major consumers taxonomic groups (cyclopoids, harpacticoids and amphipods) were assessed through Permutational multivariate analysis of variance (PERMANOVA, R-package⁴⁶ ‘vegan’) and pairwise post hoc pairwise multilevel comparisons⁴⁸.

For comparison of potential shifts in abundances of microbial metabolic pathways within groundwater samples, copepods and amphipods across rainfall periods, analysis of variance (ANOVA) was performed on the abundance data (two replicates per each group) on the predicted pathways depicting carbon fixation, carbohydrate, lipid, amino acid, methane, nitrogen and sulfur metabolisms. ANOVAs coupled with Tukey’s HSD pairwise comparisons (R-package⁴⁶ ‘stats’) were employed to inspect significant differences between bulk SIA (δ¹³C and δ¹⁵N) and essential amino acid (δ¹³C_Phe, δ¹³C_Arg, δ¹³C_Val and δ¹³C_Val-Phe) data from the stygofaunal taxa within the different rainfall conditions (LR and HR). PERMANOVAs (R-package⁴⁶ ‘vegan’) were also performed to investigate the potential clustering trends within the stygofaunal taxa across rainfall periods from the combination of radiocarbon (Δ¹⁴C) and carbon SIA (δ¹³C) isotopic fingerprints.

SIMM (Stable Isotope Mixing Models, R-package⁴⁶ ‘simmr’) were used to estimate dietary proportions of copepods and amphipods within a Bayesian framework. Due to the lack of species-specific trophic discrimination factors for stygofauna, we employed the widely accepted values of 3.4 ± 2‰ for nitrogen and 0.5 ± 1‰ for carbon⁴⁹. 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.

Source: Ecology - nature.com