Refining trophic dynamics through multi-factor Bayesian mixing models: A case study of subterranean beetles

Food web dynamics are vital in shaping the functional ecology of ecosystems. However, trophic ecology is still in its infancy in groundwater ecosystems due to the cryptic nature of these environments. To unravel trophic interactions between subterranean biota, we applied an interdisciplinary Bayesian mixing model design (multi-factor BMM) based on the integration of faunal C and N bulk tissue stable isotope data (δ13C and δ15N) with radiocarbon data (Δ14C), and prior information from metagenomic analyses. We further compared outcomes from multi-factor BMM with a conventional isotope double proxy mixing model (SIA BMM), triple proxy (δ13C, δ15N and Δ14C, multi-proxy BMM) and double proxy combined with DNA prior information (SIA+DNA BMM) designs. Three species of subterranean beetles (Paroster macrosturtensis, Paroster mesosturtensis and Paroster microsturtensis) and their main prey items Chiltoniidae amphipods (AM1: Scutachiltonia axfordi and AM2: Yilgarniella sturtensis), cyclopoids and harpacticoids from a calcrete in Western Australia were targeted. Diet estimations from stable isotope only models indicated homogeneous patterns with modest preferences for amphipods as prey items. Multi-proxy BMM suggested increased - and species-specific - predatory pressures on amphipods coupled with high rates of scavenging/predation on sister species. SIA+DNA BMM showed marked preferences for amphipods AM1 and AM2 and reduced interspecific scavenging/predation on Paroster species. Multi-factorial BMM revealed the most precise estimations (lower overall SD and very marginal beetles’ interspecific interactions), indicating consistent preferences for amphipods AM1 in all the beetles’ diets. Incorporation of genetic priors allowed crucial refining of the feeding preferences, while integration of more expensive radiocarbon data as a third proxy (when combined with genetic data) produced more precise outcomes but close dietary reconstruction to that from SIA+DNA BMM. Further multidisciplinary modelling from other groundwater environments will help elucidate the potential behind these designs and bring light to the feeding ecology of one the most vital ecosystems worldwide.

Methods
Supporting information Sacco et al., 2020 Ecology and Evolution

Extraction of DNA: a total of 15 individuals for each of the three Paroster species (P. macrosturtensis (B), P. mesosturtensis (M), P. microsturtensis (S)) were sorted into triplicate stygobitic pools (i.e. 5 individuals per pool) for DNA extraction. Prior to DNA extraction, stygobitic animals (5 individuals per pool; n=18) were placed in a petri dish containing ultrapure water and UV sterilised for 10 minutes to eliminate any potential eDNA that may be contained on the exoskeleton. Immediately post-UV treatment, the animals were placed into tissue lysis tubes with 180 µL ATL and 20 µL Proteinase K and homogenised using Minilys® tissue homogeniser (ThermoFisher Scientific, Australia) at high speed for 30 seconds. Lysis tubes, inclusive of one laboratory control, were incubated at 56°C using an agitating heat block (Eppendorf ThermoStat™ C, VWR, Australia) for 5 hours. Following incubation, DNA extraction was carried out using DNeasy Blood and Tissue Kit (Qiagen; Venlo, Netherlands eluted off the silica column in 50 µL AE buffer.

PCR: The quality and quantity of DNA extracted from each stygobitic pool was measured using quantitative PCR (qPCR), insect 16S and COI genes. PCR reactions to assess the quality and quantity of the DNA target of interest were carried out via qPCR (Applied Biosystems [ABI], USA) in 25 µL reaction volumes consisting of 2 mM MgCl2 (Fisher Biotec, Australia), 1 x 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 (MZArtF and MZArtR Zeale, Butlin, & Jones 2011; Ins16S_1shortF and Ins16S_1shortR Clarke, Soubrier, & Cooper 2014), 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 minutes, followed by 40 cycles of 95°C for 30 seconds, 52°C for 30 seconds, 72°C for 30 seconds, and a final extension at 72°C for 10 minutes.

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-8bp multiplex identifier tag (MID-tag) with the 16S and COI primer set. Independent MID-tag qPCR for each stygobitic pool were carried out in 25 µL reactions containing 1 x PCR Gold Buffer, 2.5 mM MgCl2, 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-600bp) 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 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. Each stygobitic pool for the two genes amplified were searched using BLASTn (Altschul, Gish, & Lipman 1990), against the NCBI GenBank nucleotide database to enable taxonomic identification. This was automated in the internet-based bioinformatics workflow environment through Pawsey Supercomputing (REF). The BLAST results obtained were imported into MEtaGenome Analyzer v4 (MEGAN; Huson, Auch, & Schuster 2007), where they were taxonomically assigned using the LCA-assignment algorithm (min. bit score = 65.0, top percentage = 5%, min. support = 5).

Additional references

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., & Lipman, D. J. (1990). Basic local alignment search tool. Journal of molecular biology, 215(3), 403-410.

Clarke, L. J., Soubrier, J., Weyrich, L. S., & Cooper, A. (2014). Environmental metabarcodes for insects: in silico PCR reveals potential for taxonomic bias. Molecular ecology resources, 14(6), 1160-1170.

Huson, D. H., Auch, A. F., Qi, J., & Schuster, S. C. (2007). MEGAN analysis of metagenomic data. Genome research, 17(3), 377-386.

Zeale, M. R., Butlin, R. K., Barker, G. L., Lees, D. C., & Jones, G. (2011). Taxon‐specific PCR for DNA barcoding arthropod prey in bat faeces. Molecular ecology resources, 11(2), 236-244.