Overview

We dissected digestives tracts from largemouth bass and channel catfish that were sacrificed at regular intervals post ingestion of a meal to determine how different covariates affect digestion rates. We define digestion rates as the reduction in prey item DNA detected in predator gut contents through time due to the processes of degradation, digestion, and evacuation. We used two different molecular methods to quantify the amount of DNA within the digestive tracts at each time interval and then used generalized linear models to quantify effects of multiple covariates (i.e., species, ration size, and temperature) on digestion rates.

Controlled Experiment

We conducted experiments on wild-caught and hatchery-reared largemouth bass (n fish=382) and hatchery-reared channel catfish (n fish=220) to investigate how temperature, predator species, and ration size influence digestion rates (Table 1). Hatchery largemouth bass (n = 202) and channel catfish (n = 220) were obtained from a private aquaculture facility in Galt, CA (The Fishery Inc.). The wild largemouth bass (n = 225) were captured via electroshocking in the Delta. All fish were transported to the University of California – Davis Center for Aquatic Biology and Aquaculture (UC – Davis CABA) for gastric evacuation experiments. Predator fish were transported at 20˚C, then held at 18.5˚C for three days to initiate acclimation and reduce the stress from transport. After the initial three-day acclimation period, fish were placed in their respective 3 m flow through temperature treatment tanks (15.5 ˚C and 18.5 ˚C) and cages for acclimation. The treatment temperatures were selected based on the range of temperatures observed in the Delta during the outmigration period of fall and late-fall Chinook (March – June), which are the most abundant ecotypes. We held 55 predators within each treatment tank, for a total of 110 predators per trial. Acclimation for two weeks was deemed necessary to 1) ensure previously digested material was eliminated and, 2) reduce overall stress and allow digestion rates to return to  normal levels (Brandl et al., 2016). Ration sizes of one and three smolts were chosen based on the typical amount of prey found within a piscivore’s stomach (Nobriga & Feyrer, 2007). Within each treatment tank, fish were placed into 4 cages to organize the predators based on time of ingestion throughout the feeding trials. We started the trials with 5 predators for each combination of temperature and ration (e.g., 15 ˚C and 1 smolt), thus, ideally we would have a total of 200 samples throughout 10 dissection intervals (see below). Each tankcontained 5 spare predators in the case of mortality or regurgitation. 

Prey items for the diet study were fall run Chinook salmon smolts received from state run hatchery in  elements, California (Mokelumne hatchery). Smolts were reared in 1.2 m flow through tanks set to 12˚ C and fed fish pellets through automated feeders. Smolts were sacrificed by a single strike to the top of the skull, wet weighed (g), fin clipped for DNA identification, and force fed to predators. For the feeding trial with hatchery-reared largemouth bass, smolts were divided into smaller portions due to smaller stomach sizes of the hatchery-reared predators. Smolts were cut near the posterior portion of the gill operculum and fish were fed the ⅓ portion that included the head while the remainder of the smolt was discarded. After a two-week acclimation period, predators were anesthetized in a buffered solution of MS-222, tagged for identification (Floy tags, Seattle, Washington), measured to nearest fork length (mm), and force fed their respective smolt ration. Feeding was conducted by placing euthanized smolts down the pharynx of immobilized predators using forceps. Once fed, predators were promptly returned to their respective treatment tanks and cages for recovery. During the first 30 minutes of recovery, predators were observed for regurgitation. If regurgitation occurred, the predator sample was removed and a new predator from the acclimated spares was fed and added to the tank. If regurgitation was suspected beyond the 30-minute observation period, a spare predator was force-fed and added to that trial as a replacement. Throughout all experimental treatments, a total of 10 regurgitations and 8 premature mortalities were observed (Supplemental Table 1). Regurgitations occurred across all species, temperatures, and feed ration sizes. Although regurgitation occurred across all trials, more occurred within warm water treatments (7 in 18.5°C vs 3 in 15.5°C). We also note that premature mortalities were observed only on wild largemouth bass, and not on hatchery subjects of either species. We hypothesize that the stress associated with electroshocking and transport affected the survival of wild largemouth bass.

At regular intervals post-ingestion, a subset of five predators from each treatment were euthanized by a single strike to the top of the skull and the digestive tracts were removed. The entire digestive tract was considered a single sample and consisted of the trachea, the pyloric cecum, stomach, pyloric sphincter, intestinal tract, and anal cavity. To remove these samples, a sterile surgical scalpel was used to make an incision from the predator’s anal cavity along the ventral side of their body, anterior to the gill plate. A second sterile scalpel was then used to make incisions internally until the entire digestive tract could be removed. Digestive tract samples were injected with a 3ml solution of 100% non-denatured ethanol (EtOH) to halt enzymatic digestion and stored in 200 ml conical vials. Conical vials were filled with 100% non-denatured EtOH ensuring that the sample to liquid ratio was approximately 2:3. Sample vials were refreshed with 100% non-denatured EtOH 24 hours post dissection and wrapped in parafilm to prevent evaporation. Samples were then stored at room temperature until DNA extraction began. To investigate the persistence of ingested prey over time, study subjects were sacrificed, and gut contents were sampled at regular intervals after the feeding treatments, with the earliest dissections commencing 6hrs post ingestion followed by regularly spaced dissections at 12h, 24h, 36h, 48h, 60h, 72h, 84h, 96h and a final dissection period at 120h. A total of 50 predator samples per species were targeted per temperature per 231   ration size.

Visual diet analysis

After isolating the diet samples, we conducted a visual assessment to quantify the relative amount of prey material within the digestive tract. The two goals of this assessment were: (1) to determine if the quantity of material in the digestive tract was correlated with prey DNA quantity, and (2) to determine if the quality of digestive tract material (i.e., the ratio of tissue vs. stool) affected the molecular analyses. To conduct the visual assessment, we first removed the prey contents from the entire digestive tract (stomach and intestines) and preserved them in ethanol. While removing the prey contents, we noted whether stool was present or absent and whether it was a low, medium, or high amount relative to other samples (rank from 0 for absent to 3 for high). We used the same scoring system for undigested tissue. Each sample could therefore have a score from 0 (no stool or tissue present) to 6 (high stool and high tissue). Stool was defined as relatively large and generally dark in color, while undigested tissue was lighter in color and in chunks of variable size and shape. For undigested tissue, we also noted whether

we could identify the prey item’s taxa based on external characteristics. After visual assessment, we prepared each stomach for DNA extraction. First, we pipetted as much stool as possible into a 1.5 ml tube. If no stool was visually present, we still pipetted approximately 700 µl of ethanol from the bottom of the vial used to preserve the stomach into the tube for extraction. To sample chunks of partially digested tissue, a small piece was taken from each visible chunk to ensure thorough sampling. We recognize that our sampling approach is not fully quantitative but it represented a logistically feasible approach that should be semi-quantitative. DNA was extracted from this mixture of stool and tissue (see below).

Molecular Methods

Before extraction, stomach and intestinal contents were centrifuged for 3 minutes at 5,724 x g and any ethanol was removed by pipette and discarded. Excess ethanol was then allowed to evaporate overnight. Extractions were conducted with DNA stool kits (Macherev Nagel Necleospin 96) modified by replacing bead-induced lysis with enzymatic lysis, using a per-sample volume of 25 µl proteinase-k with 850 µl lysis buffer ST1. Samples were incubated overnight at 56°C, and the subsequent extraction followed the standard manufacturer protocol. One negative control was included on each 96-well extraction plate. It is important to note that the samples were extracted and PCR was conducted in a laboratory space where Chinook salmon tissue is extracted and genotyped for nuclear SNPs (PCR and extraction rooms are kept separate). We are confident that contamination did not influence our findings as it is clear that our molecular results reflected known biological processes (i.e. DNA detectability decreases as digestion progresses). However, our laboratory setup would not be appropriate for a study where stomachs from wild populations were analyzed because contamination could significantly impact results.   

All extracted stomach samples were analyzed with qPCR using an assay designed to quantify Chinook salmon DNA (Brandl et al., 2015). To prepare a standard curve, extracted genomic DNA from Chinook salmon available from Shi et al. (in press) was quantified using an assay kit (Qubit 1X HS dsDNA kit), visualized on a 2% agarose gel to determine average fragment size and calculate copy number, and then serially diluted. Each 96-well plate included seven standards in triplicate and three no-template-controls. Reactions were conducted in 10 µl volumes each containing 1x 5 µl master mix (Applied Biosystems TaqMan Universal PCR Master Mix), 0.9 µM concentration of each primer, 0.7 µM concentration of probe, and 3 µL template DNA. Thermocycling was performed using a real-time PCR system (QuantStudio 12K Flex) with the following profile: 10 min at 95 °C, and then 40 cycles of 15 seconds at 95 °C for 15 s and 60 °C for 1 minute. All samples were run in triplicate and the mean copy number was calculated for each sample. Amplification was considered successful if the Cq > 0 and amplification above baseline was registered for an unknown sample using the PCR system software (QuantStudio). If some replicates did not show amplification for a given sample the mean was calculated using the replicates that did amplify.

We also conducted metabarcoding on all diet samples using the MiFish 12S primer set (Miya et al., 2015). No PCR replicates were conducted for this analysis. Libraries were prepared with a two-step PCR similar to the GT-seq protocol outlined in Campbell et al. (2015). Initial PCR reactions were performed in 12 µl volumes using 2 µl of template DNA and 10 µl of PCR master mix. The PCR master mix consisted of per reaction volumes of: 6 µl master mix (Qiagen Multiplex PCR Plus), 3.4 µl sterile water and 0.3 µl each of the MiFish forward and reverse primer at 10 µM. Thermal cycling was performed as follows: 95°C for 5 minutes, followed by 35 cycles of 95°C for 30 seconds, 65 °C for 15 seconds, and 72°C for 15 seconds, an extension at 72°C for 5 minutes, and a 8°C hold.  PCR products were indexed in a barcoding PCR as dictated by Campbell et al. (2015), modified by reducing working i05 primer concentrations to 0.5 µM and working i07 primer concentrations to 1 µM. Following this, each 96-well plate was subsequently pooled without normalizing the amount of DNA per sample to maintain a relationship with the initial concentration of DNA in the extract. A double-sided bead size selection was performed using paramagnetic beads (Beckman Coulter AMPure XP beads), using ratios of beads to library of 0.5x to remove non-target larger fragments and then 0.7x to remove dimer and retain the desired amplicon. Libraries were sequenced with a next-generation sequencing platform (MiSeq Illumina) using a single v2 300-cycle kit with 2 x 150 bp  aired-end (PE) chemistry with 15% PhiX to compensate for the low diversity library. The lane included  82 samples loaded at equal concentrations, 163 of which were not part of the current study, and 16 of which were no-template-controls (PCR negatives). Sequencing reads were assigned to individual samples using the sequencing platform analysis software (Illumina MiSeq Analysis Software). Primers were trimmed from forward and reverse reads using open-source software (Cutadapt) and reads without primer sequences were discarded (Martin, 2011). Primer-trimmed fastq files were imported into open-source software (DADA2; Callahan et al., 2016) in R for quality filtering, chimera removal, and identifying amplicon sequence variants (ASV) using the pseudo-pool option for the core dada2 algorithm (additional parameters included truncLen = c(110,100), maxN=0, maxEE=c(2,4) and truncQ=2). Forward and reverse reads were merged and then off-target sequences that likely amplified non-fish amplicons were removed using a length filter of 166-172 bp based on the expected amplicon size for MiFish (with primers removed). Output files included a table of the number of reads per ASV per sample (ASV table) and the sequence identity for each ASV, which was used as input for a nucleotide database search (BLASTn; Camacho et al., 2009) using a local copy of the National Center for Biotechnology Information (NCBI) nucleotide database and a minimum 96 percent sequence identity (for the entire 166-172 bp amplicon). The software program taxonkit (Shen & Ren, 2021) was used to match taxonomic lineage to the BLASTn results. This taxonomy information was imported into R and custom scripts were used to assign ASV sequences to the appropriate taxonomic level. When multiple species matched the same ASV within 2% sequence identity of one another, taxonomic levels (i.e., genus, family, order) were iteratively increased until matches were unambiguous. Species-level assignments were only accepted at 98% identity, whereas matches to all other taxonomic levels were retained at 96% sequence identity. This taxonomic information was combined with samples and read count data. The minimum threshold for retaining metabarcoding data was two reads per taxon per sample. We recognize that this value is lower than many empirical studies but we felt it was appropriate because we wanted as much data available to model and because this was a controlled experiment rather than an empirical study with the goal of detecting prey items in wild populations.

Data Analysis

Correspondence of visual assessments and molecular gut content analysis We hypothesized that the quality of the prey DNA would degrade over time through the process of digestion. To test this hypothesis, we used Pearson correlations to determine if the quantity of diet contents (stool + prey tissue scores) was correlated with smolt prey item read counts (metabarcoding) or the log transformed copy number (qPCR) from molecular analyses. We analyzed each species separately, but we pooled gut samples across temperature and ration treatments to ensure we had a sufficient sample size at each dissection time interval post prey ingestion.

Regression modeling of digestion rates

We used a generalized linear model (GLM) with a Tweedie distribution to determine what factors were most correlated with the abundance of prey item DNA in sampled guts. Models were fit separately for the two molecular methods, with read count the response variable for metabarcoding diet analysis and copy number the response variable for qPCR analyses. The responses of both methods contained a substantial number of zeros values as prey DNA was digested away or evacuated. Furthermore, both qPCR and metabarcoding were implemented using minimum genetic material thresholds that need be exceeded to generate positive DNA detection (> 2 reads for metabarcoding, Cq > 0 for qPCR; see above). Thus, we used GLMs with a Tweedie distribution, which is a special case of the exponential distribution that is very flexible and can account for data that contains a large number of zeros (Tweedie, 1984). Tweedie GLMs also have an advantage over zero-inflated models because the zeros are handled uniformly, which improves the interpretability of model coefficients (Shono, 2008). The package “glmmTMB” (Brooks et al., 2017) was used within R (R, Development Core Team, 2021) to fit the GLM. Due to the non-normality of our model, we used quantile residuals simulated using the “DHARMa” package (Hartig, 2022) to assess model diagnostics. 

We selected four factors (time post-injection, species, temperature, and meal size) to include in our model based on a priori hypotheses regarding their relationship with digestion rates. The first factor was time post-ingestion, which is expected to reduce DNA quantity due to the physical and chemical digestion of prey tissue (Moran et al., 2016). The second factor was species, because there are numerous behavioral and physiological differences betweenlargemouth bass and channel catfish. These two species have substantially different feeding ecologies, but similar temperature preferences. Both largemouth bass and channel catfish prefer relatively warm temperatures in the 25-30 °C range, but can tolerate a large range of temperatures (Venables et al., 1978; Wellborn, 1990). Largemouth bass are visual predators and become piscivorous as early as 5 cm in length, leading to a predominantly piscivorous life history (Davis & Lock, 2007). In contrast, channel catfish are generally classified as omnivores, feeding on fish but also plants, invertebrates, amphibians, and other organisms (Braun & Phelps, 2016; Hill et al., 1995). Carnivores typically have shorter intestines and faster digestion rates, whereas herbivores and omnivores have longer digestive tracts and slower digestion rates (reviewed in Volkoff & Rønnestad, 2020). The third factor in our model was temperature, which can affect digestion rates differently depending on gut morphology (Volkoff & Rønnestad, 2020). Therefore, as water temperatures change throughout a season, digestion rates in largemouth bass and channel catfish may not respond in the same way. In our experiment, water temperature was designated as a binary factor based on the two average temperatures of the treatment tanks: 15.5˚C and 18.5˚C. Higher tank temperature is expected to increase the metabolic rate of fish and subsequently increase evacuation rate (Brett & Groves, 1979). In contrast, cooler temperatures are expected to slow digestion rates, potentially prolonging DNA detection. The last factor in our model was meal size, which could influence digestion rates due to longer processing times for larger prey items (Beamish, 1972). The binary rations of 1 and 3 smolts were transformed to a continuous factor called predator-to-prey ratio (PPR) because we hypothesized that the relative size of the prey was more informative to the digestion rate than the ration size (Supplemental Figure 1). PPR has been used frequently in food web and community sized predation studies to understand metabolic constraints of predators (Woodward & Warren, 2007). We hypothesize that we would observe an inverse relationship between PPR and digestion rate, whereby a larger PPR would decrease gastric processing time. To calculate PPR, smolts were wet weighed to nearest 0.1 g and predator weights were estimated based on length-weight regressions. Predators fork lengths were measured to nearest 1 mm and weight-at-length regressions from the literature were used to estimate individual predator weights (Henson, 1991; Keenan et al., 2011). We used the following weight-at-length regression for largemouth bass:

                                                                              ?(?) = (??)?

 where weight in grams (w) at length in mm (l). An alternative formulation of the model was used for channel catfish:

                                                                          ?(?) = (L)lb 

where L is equal to the standard length of a catfish at 1 kg (Keenan et al., 2011). We used this alternative formulation for Channel catfish because they are known to have a lean, longer form at higher weights compared to largemouth bass. No morphological differences were observed between the two sources of largemouth bass, so the same length-weight regression was used on hatchery and wild largemouth bass. We then divided the calculated predator mass by their respectively fed prey mass to calculate the PPR. Higher PPR means that the predator was large relative to the size of the meal. In other words, a predator of a given size that ate a small meal would have a higher PPR than the same predator that ate a large meal.

After calculating this ratio, we determined that the PPR was correlated with predator weights and that 14 catfish had much higher PPRs (>300) than all other fish in the study (Supplemental Figure 1). Furthermore, these large catfish were all in the single fish ration treatment. To better compare evacuation rates between species based on the experimental treatments, we excluded these 14 large catfish from our modeling analysis.

We assessed support for different model structures in our regressions of feeding trial factors against molecular diet analysis assay outcomes (read counts for metabarcoding and copy number for qPCR) using Akaike’s Information Criterion (AIC; Burnham, 1998). Our primary model selection goal was to identify the top AIC-supported regression amongst a larger candidate model set made up of plausible model structures. Time post-ingestion was expected to have the largest effect on the measured quantity of DNA, and therefore was included in each candidate model. Temperature, PPR, and species were included in models as main effects and as an interaction with time post-ingestion. The interactions between time and each of the other covariates were included to determine if there was an effect of any covariate on the digestion rate. We also included an interaction term between species and PPR because we hypothesized that ration size could influence ingested prey evacuation rates differently across species due to differences in metabolism and digestive morphology. Finally, we included an interaction between species and temperature because largemouth bass and channel catfish have different physiologies. Thus, their prey digestion rates may not change proportionally to changes in temperature. Interaction terms were only included if the variable was included in the model as a main effect. We conducted model selection using the “MuMIn” package (Barton, 2022).

Genetic half-life of prey DNA

While our GLM regressions on molecular diet analysis read counts (metabarcoding) or copy number (qPCR) characterize the influence of feeding trial factors on the persistence of detectable prey item DNA, we also calculated prey DNA genetic half-life values to provide a convenient single metric to facilitate DNA loss comparisons across samples. Detectability half-life provides a measure reflective of the rate that detectable prey item DNA is lost from a predator gut as prey rations are digested, degraded, and evacuated from the gut tract. Wedefine detectability half-life as the duration it takes for smolt DNA in gut content samples to reduce to half of its initial value. To calculate the genetic half-life (HL) of prey item DNA in predator gut content samples, we used the following equation:

                                                                             HL =    ??

                                                                                     ?????,?

                                                                                  ???0.5 ????

Where t is amount of time that passed between interval (i) and interval (j) and ?DNA is the difference in DNA quantity between interval (i) and interval (j). DNA quantities at intervals i and j were predicted from the respective top AIC-supported GLM regression models fit to read count (metabarcoding) or copy number (qPCR) molecular diet analysis data described above (see:‘Regression modeling of prey DNA loss’). Because they were estimated based on the model, we will refer to them as ‘GLM-estimated half-life’ values. Prey item DNA half-life values were summarized across all combinations of feeding trial factors (species, temperature, and PPR). We used parametric bootstrapping (1000 samples) to estimate the 95% confidence interval for our half-life estimates based on the mean and standard deviation of the predicted DNA quantity. To further assess the relative strength of influence of different feeding trial factors on digestion rates we calculated a half-life ratio. The half-life ratio (HLrat) was calculated using the equation:

                                                                                          ??k

                                                                              ????? = ??l

Where HLk is the half-life value for one factor at level k (e.g., species = CCF) and HLl is the half-life value for the same factor at level l (e.g., species = LMB) when the other two covariates are equal (e.g., temperature = 15.5 & PPR = low). This was repeated for every combination of factors.