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

Abstract 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 (SIA BMM) 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 modeling 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.


| INTRODUC TI ON
Trophic dynamics provide vital information about ecological functioning (Lindeman, 1942;Polis & Winemiller, 2013;Start, 2018). Food web interactions shape ecological niche occupations and frame community dynamics (de Ruiter, Wolters, Moore, & Winemiller, 2005). The functionality of the trophic web is ultimately defined by intra-and interspecific interactions which shape biochemical patterns and energy flows within the ecosystems (Begon, Townsend, & Harper, 2006).
Both qualitative (Paine, 1980) and quantitative (Banasek-Richter, Cattin, & Bersier, 2004) approaches have been applied in several ecosystems, the latter being more accurate but more challenging than the former (Kadoya, Osada, & Takimoto, 2012). Over the last four decades, isotope mixing models, such as IsoSource (Phillips & Gregg, 2001) or Bayesian mixing models (BMM, Parnell et al., 2013), have been increasingly used for quantitative reconstruction of dietary preferences. Both techniques aim to quantify unknown mixing contributions via measurement of the isotopic signals in consumers and food sources (Post, 2002).
Dietary proxies based on bulk carbon (δ 13 C) and nitrogen (δ 15 N) stable isotope analysis (SIA) are a powerful tool for studying trophic preferences and food web interactions (Fry, 2006 and references therein). Concurrently, radiocarbon ( 14 C) forms a key tracer in untangling carbon incorporation and trophic pathways (Larsen, Yokoyama, & Fernandes, 2018). Metagenomics data, integrated with consumer and source abundances, provide semi-quantitative prior information on dietary preferences that can refine statistical modeling (Chiaradia, Forero, McInnes, & Ramírez, 2014). BMM (i.e., MixSiar (Stock & Semmens, 2016)) allows integration of data from different disciplines such as biochemistry, genetics, and ecology, but the majority of trophic studies focus on stable isotopic frameworks.
Isotopic analysis (δ 13 C and δ 15 N SIA) of the three diving beetles revealed that the predatory pressures on both amphipods and copepods were also coupled with marginal interspecific predatory pressures on Paroster species (Saccò, Blyth, Humphreys, et al., 2019). However, the analysis of trophic interactions through such conventional approaches faces major challenges (Boecklen, Yarnes, Cook, & James, 2011), stressing the need for cost-efficient model designs that allow the combination of data from multiple disciplines (Saccò, Blyth, Bateman, et al., 2019).
Here, we test whether the multi-factor design of FRUITS models enables refinement of dietary preferences in the three species of subterranean aquatic beetles along with their food sources. The work aims to (a) evaluate the use of multi-discipline and/or isotope only models in subterranean ecosystems and (b) provide recommendations on the use of isotopic techniques in groundwater ecology.

| Fieldwork, sample preparation, and faunal trophic ecology
Sampling occurred at a calcrete aquifer on Sturt Meadows pastoral station, in the Yilgarn, Western Australia. Stygofauna were sampled by haul nets (100 µm mesh size) from boreholes during two sampling campaigns carried out in July and November 2017. For further details about the sampling design, see Saccò et al. (2020). Specimens were sorted under a stereomicroscope to species level with reference to specific taxonomic keys (King, Bradford, Austin, Humphreys, & Cooper, 2012;Watts & Humphreys, 2009). All individuals from a single taxon were combined into one pool and washed with MilliQ water to remove external contaminants. Samples were oven dried

K E Y W O R D S
Bayesian mixing models, food webs, groundwater, metagenomics, radiocarbon, stygofauna at 60°C overnight, crushed to a fine powder, and stored at -20°C until analysis.
More details about the taxa morphology are provided in the Note S1.
Investigations on trophic habits based on genetic information (metabarcoding analysis) from previous studies at SM calcrete suggested that all three Paroster species feed on amphipod AM1 (ranging from 68% (B) to 28% (S) of their diet) more than the other groups and avoid intraspecific cannibalism. Specifically, while the diet of beetles B was dominated by amphipods (90% overall), beetles M and S preferred harpacticoids over amphipods AM2 (Bradford et al., 2014) (Table 1). Multi-primer metabarcoding analyses (Ins16S and MZartCOI) on the three species indicated that occasional reciprocal scavenging/predation on sister species follows sister species-specific patterns (Hyde, 2018) (Table S1).
Estimation of diet proportions of B, M, and S from isotopic analyses confirmed these trends, with amphipod AM1 being the preferred prey for all the three predator species (B: 25%; M: 27.4%; S: 25.4%) and copepods accounting for the 30% beetles' diets (Saccò, Blyth, Humphreys, et al., 2019).

| Bulk SIA
C and N bulk stable isotopic analyses on homogenized stygofaunal samples (1.28 mg, 0.08-0.14 mg, and 0.63-2.79 mg per samples, respectively, see Saccò, Blyth, Humphreys, et al. (2019) for further details) were performed at the Australian Nuclear Science and Technology Organisation (ANSTO). Samples were loaded into tin capsules and analyzed by a continuous flow isotope ratio mass spectrometer (CF-IRMS), model Delta V Plus (Thermo Scientific Corporation), interfaced with an elemental analyser (Thermo Fisher Flash 2000 HT EA, Thermo Electron Corporation) following Mazumder, Saintilan, Wen, Kobayashi, & Rogers, (2017). Carbon and nitrogen isotopic values are reported in per mil (‰) according to the standard delta (δ) notation, relative to the Vienna Peedee Belemnite (VPDB) and to atmospheric nitrogen (AIR), respectively. Results have an analytical precision of ±0.30‰. Results on the % of C and %N from bulk tissue were also obtained through the elemental analyser.

| Radiocarbon
Stygofaunal samples (~1 mg per sample for beetles (B, M and S) and amphipods (AM1 and AM2) and ~0.5 mg for copepods) were treated with dilute HCl (1 M) for 2 hr to remove carbonate contamination.
Due to sample size constraints, cyclopoids and harpacticoids were combined in one sample, and therefore, a unique radiocarbon value (the first ever recorded for groundwater copepods) for both groups was obtained. The pre-treated samples were combusted to CO 2 and converted to graphite following Hua et al. (2001). 14 C content of the sample graphite was determined using the accelerator mass spectrometry (AMS) STAR Facility at ANSTO (Sydney, Australia; Fink et al., 2004). Radiocarbon results are reported in Δ 14 C value in per mil (‰) relative to the absolute radiocarbon standard activity, and age was also assessed (with present being 1950 AD) (Stuiver & Polach, 1977),

| Statistical analysis
Relative contributions from dietary items were estimated using the software FRUITS version 2.1.1 beta (Fernandes et al., 2014). FRUITS

F I G U R E 1 Photographs illustrating specimens belonging to the species (B) Paroster macrosturtensis (a), (M) Paroster mesosturtensis
allows quantification of the dietary proportions of food sources (defined as "sources") among consumers (defined as "targets") via isotopic quantitative signals (defined as "proxies"). The model incorporates variance associated with isotopic measurements from sources and targets, trophic discrimination offsets, and allows incorporation of prior information to refine the analysis (defined as "priors"). FRUITS models generate a BUGS (Bayesian inference Using Gibbs Sampling) coding that is then transferred to OpenBUGS package, a software commonly used for Bayesian probability modeling (Lunn, Thomas, Best, & Spiegelhalter, 2000). Markov chain Monte Carlo (MCMC) simulations allow generation of posterior distributions associated with credible intervals (Gilks, Richardson, & Spiegelhalter, 1996). BMM was applied to test the potential changes in dietary preferences. Figure 2 depicts the statistical design used for the present study.
The multi-factor BMM was run using δ 13 C and δ 15 N values from bulk tissue analysis. Specific trophic discrimination factors are not yet available for stygofauna, so in all bulk tissues we used the widely accepted discrimination values of 3.46 ± 2‰ for nitrogen and 0.5 ± 1‰ for carbon (Zanden & Rasmussen, 2001). The third proxy was radiocarbon data (Δ 14 C). Coupled with stable isotope proxies, this has the potential to provide more discriminatory carbon fingerprints between sources and consumers (Larsen et al., 2018). Δ 14 C is free from isotopic fractionation due to the internal correction by the δ 13 C of −25‰ (Ishikawa, Hyodo, & Tayasu, 2013). This "triple proxy approach" (δ 13 C, δ 15 N, and Δ 14 C) was coupled with inputs from prior information (both from sources and targets) from metagenomics analyses to create a multi-factorial design (hereafter defined as "multi-factor BMM").
Despite the differences in research targets between isotope ecology (focused on "who assimilates whom") and metagenomics (focused on TA B L E 1 Trophic behaviors (predation and scavenging) based on prior metagenomics data on Paroster macrosturtensis (B), Paroster mesosturtensis (M), and Paroster microsturtensis (S)

| Multi-proxy SIA BMM
Diet estimations of beetles B were markedly dominated by amphipods AM1, which composed almost two thirds of the overall con-  Figure 3a.2). The diet of beetles S was dominated by amphipods AM1 (62.60 ± 12.68%) and AM2 (11.07 ± 8.85%), with sister species scavenging/predation sitting at the same secondary level as copepods' predation (Figure 3a.3).

TA B L E 2
Carbon and nitrogen isotopic ratios for δ 13 C and δ 15 N SIA (together with the % of C and N) and Δ 14 C values (in ‰)

| Multi-factor mixing models in calcretes
Our  (Bradford et al., 2014) and ecological (Hyde, 2018) forces are expected to drive very marginal interspecies predation among top predators at SM calcrete. As a result, we argue here that the outcomes from both SIA and multi-proxy BMM provide imperfect beetles' diet estimations at SM calcrete and advocate for the integration of further data for more reliable modeling.
Compared to previous designs, the multi-factorial model (multi-factor BMM) showed reduced levels of uncertainty (SD values always below 10%, Table S2) and suggested that subterranean beetles exert preferential predatory pressures on AM1 coupled with extremely marginal (always below 5% combining two groups) interspecific interactions.
Our outcomes from multi-factor modeling, together with the results from SIA + DNA BMM, indicated that subterranean beetles at Sturt Meadows lack trophic niche partitioning, as already suggested by Bradford et al. (2014). While seemingly counterintuitive and in contrast to the classic subterranean ecology paradigm of opportunistic feeding traits, our results coincide with the conclusions drawn by Francois et al. (2016). This indicated that widely reported low metabolic rates and resource-gathering abilities might play a role in releasing the constraints on trophic specialization underground. Resulting from these eco-evolutionary dynamics, groundwater fauna is suggested to display feeding habits focused on more ubiquitous resources (i.e., sedimentary biofilm, prey items, etc.) rather than being driven by selective forces toward generalist strategies. At SM calcrete, species-specific ethological (i.e., group feeding) and physiological features (high efficiency in metabolic activation processes in the smaller species such as M and S, see 3) for using bulk tissue δ 13 C and δ 15 N combined with prior metabarcoding data (SIA + DNA BMM). Boxes and whiskers indicate 68% and 95% credible intervals, respectively. Horizontal continuous lines indicate the estimated mean, and dashed lines refer to the median. Refer to Table S2 for the specific contribution mean values beetles. However, further work is needed to determine which model is most accurate. Indeed, prior information in the posterior distribution is likely to have played a key role in shaping our patterns, confirming that estimation of diet benefits from genetic data on potential prey (Chiaradia et al., 2014). Given that the majority of information about groundwater fauna comes from genetic investigations (i.e., metabarcoding), BMM allows novel coupling between isotopic (quantitative) and molecular (qualitative or semi-quantitative) data and has the potential to bring light to the mechanisms sustaining biodiversity in of one of the largest and most understudied ecosystems in the world. This combination of methodologies remains under-exploited (Majdi et al., 2018), with the present study still limited to one arid zone calcrete system. Investigations involving different groundwater environments (i.e. alluvial aquifers, karst, etc.) will help elucidate the potential behind the integration of data from different disciplines into isotopic ecology studies on subterranean fauna.

| Defining the (realistic) current isotopic design in groundwater studies
The study of stygofaunal foraging ecology through isotopic techniques is gaining prominence as an analytical approach to dig into groundwater energy flows and trophic niche interactions (Saccò, Blyth, Bateman, et al., 2019). However, the current technical and analytical advances seen in the broader field of isotopic ecology are frequently coupled with increased price. As a result, a balance between cost and precision of outcome must be achieved.
To date, the vast majority of groundwater isotopic food web studies involve conventional bulk tissue SIA (e.g., Hartland, Fenwick, & Bury, 2011;Simon, Benfield, & Macko, 2003). However, δ 13 C and δ 15 N measurement alone, despite being the cheapest analytical approach available, have been reported to be only partially accurate due to the mixing of biochemical fractionation pathways (Newsome, Fogel, Kelly, & del Rio, 2011 and references therein). Our results concur with this observation, indicating that isotopic trophic studies in groundwaters using classic SIA designs are potentially exposed to misinterpretation ( Table 3).
The incorporation of a third proxy (Δ 14 C in our study) into BMM is unexplored in groundwater feeding ecology studies, probably due to budgetary constraints ( 14 C analysis has a cost that can exceed ten times SIA) and analytical issues ( 3) for using bulk tissue δ 13 C and δ 15 N, combined with Δ 14 C and prior metabarcoding data (multi-factor BMM). Boxes and whiskers indicate 68% and 95% credible intervals, respectively. Horizontal continuous lines indicate the estimated mean, and dashed lines refer to the median. Refer to Table S2 for the specific contribution mean values beetles' interspecific interactions indicated by the previous studies carried out at SM calcrete (Bradford, 2010;Bradford et al., 2014;Hyde, 2018). Overall, this scenario suggests that the only partial increase in precision does not merit the additional cost of the triple-proxy design without additional prior information (Table 3).
Stygofauna can display cryptic feeding habits (Stoch, 1995) which are hard to investigate in mesocosm experiments, and genetically based biological characterizations can be crucial in identifying diet preferences under natural conditions (Saccò, Blyth, Bateman, et al., 2019). Once incorporated into BMM, this combination of data from independent sources enables comparison of datasets and can constrain the systemic biases of each separate technique (Chiaradia et al., 2014). Overall, we suggest that reconciliation of the trade-off between the cost (price) and benefit (precision) in groundwater food web studies can be achieved by incorporating metabarcoding data into a model with conventional SIA. The most precise outcomes are obtained by integration of triple isotope proxy and DNA data, but overall dietary reconstruction was close to that from bulk tissue isotopes coupled with genetic prior information. Therefore, SIA + DNA is recommended generally, with full multi-factorial approaches used where operational costs are not a significant constraint.
Further advances, including specific investigations on the variability of the trophic discrimination factors (McMahon & McCarthy, 2016) for stygofauna, will enhance the biochemical understanding of trophic pathways and help refine analyses (Parnell, Inger, Bearhop, & Jackson, 2010). Recent more expensive novel analytical approaches such as compound specific isotopic analyses offer to refine foraging ecology studies and overcome some of the homogenization issues in bulk tissue SIA (Chikaraishi et al., 2009;Larsen et al., 2013;Steffan et al., 2013). The combination of SIA and CSIA has recently gained prominence in the broad literature (Potapov, Tiunov, Scheu, Larsen, & Pollierer, 2019) and has been applied in the field of groundwater ecology (Saccò, Blyth, Humphreys, et al., 2019). However, while these techniques are a cornerstone in trophic studies, conventional SIA approaches are likely to be widely used in the near future due to constraints of budget and technical limitations in CSIA. Despite its averaging of biochemical fractionation pathways, bulk tissue SIA, when integrated with prior qualitative information on the feeding habits, still allows elucidation of the food web interactions. When applied to groundwater ecology studies, we believe that these designs have the potential to enable affordable and reasonably accurate interpretation of the stygofaunal foraging ecology.

ACK N OWLED G M ENTS
We wish to acknowledge the traditional custodians of the land at

DATA AVA I L A B I L I T Y S TAT E M E N T
All additional data are available in the Supporting Information and will be archived in the Dryad repository (https://doi.org/10.5061/ dryad.2z34t mpj7).