High‐throughput sequencing outperforms traditional morphological methods in Blue Catfish diet analysis and reveals novel insights into diet ecology

Abstract Blue Catfish Ictalurus furcatus are an invasive, yet economically important species in the Chesapeake Bay. However, their impact on the trophic ecology of this system is not well understood. In order to provide in‐depth analysis of predation by Blue Catfish, we identified prey items using high‐throughput DNA sequencing (HTS) of entire gastrointestinal tracts from 134 samples using two genetic markers, mitochondrial cytochrome c oxidase I (COI) and the nuclear 18S ribosomal RNA gene. We compared our HTS results to a more traditional “hybrid” approach that coupled morphological identification with DNA barcoding. The hybrid study was conducted on additional Blue Catfish samples (n = 617 stomachs) collected from the same location and season in the previous year. Taxonomic representation with HTS vastly surpassed that achieved with the hybrid methodology in Blue Catfish. Significantly, our HTS study identified several instances of at‐risk and invasive species consumption not identified using the hybrid method, supporting the hypothesis that previous studies using morphological methods may greatly underestimate consumption of critical species. Finally, we report the novel finding that Blue Catfish diet diversity inversely correlates to daily flow rates, perhaps due to higher mobility and prey‐seeking behaviors exhibited during lower flow.

with the relative abundance of this omnivorous species presents a concern for management agencies as Blue Catfish are known to prey on native species such as economically important Blue Crab Callinectes sapidus and threatened alosines (Iwanowicz et al., 2019;Schmitt et al., 2017). Additionally, Blue Catfish can outcompete native species for resources, and White Catfish Ameiurus catus populations have declined as Blue Catfish populations have increased (Chesapeake Bay Program, 2020).
Recent studies have sought to address the full impact of Blue Catfish feeding ecology in its non-native range using visual, morphometric methods (Waters et al., 2004), sometimes in combination with DNA barcoding of specific tissues (Aguilar et al., 2017;Moran et al., 2015;Schmitt et al., 2017). Such studies are difficult to use in quantifying the effect of Blue Catfish predation on ecosystems as taxonomic identification of prey items is highly dependent upon tissue degradation, thus limiting prey identification to those items only recently consumed and of sufficient mass to withstand rapid digestive processes (Rees et al., 2020;Su et al., 2018). Studies indicate that soft-tissue prey, including larval and juvenile fishes, are unrecognizably digested within a span of twenty minutes to a few hours of consumption (Bromley, 1994;Carreon-Martinez et al., 2011).
The use of high-throughput sequencing (HTS) for diet studies is becoming more widely accepted (Casey et al., 2019;Rees et al., 2020;Waraniak et al., 2019). Such methods perform better than morphological observation, especially when diet items, including at-risk or economically important species, are partially or fully digested with no assemblance of an organism for viewing, as mentioned above (Schwarz et al., 2018;Tverin et al., 2019;Waraniak et al., 2018). Researchers can describe diverse prey assemblages (Bessey et al., 2019;Sousa et al., 2016), analyze non-native species' diets in newly invaded habitats (Harms-Tuohy et al., 2016), and quantify trophic interactions (Casey et al., 2019) using HTS methods. Despite the benefits of HTS to morphological studies and the importance of Blue Catfish from both a conservation and economic viewpoint, only one study using HTS to analyze the diet of Blue Catfish in a limited numbers of samples (n = 12) has been published to date (Iwanowicz et al., 2019).
In this study, we sought to gain a holistic view of diet behavior and of the potential impact of invasive Blue Catfish predation in a Chesapeake Bay tidal river ecosystem using high-throughput sequencing methods. To accomplish this objective, we conducted HTS of 134 Blue Catfish representing multiple life stages that were collected in the fall of 2016 from the Pamunkey River, Virginia.
High-throughput sequencing was employed on material collected from the whole gut, stomach to anus, using 18S and COI, genes commonly employed in diet analyses (Leray & Knowlton, 2015;Zhan et al., 2014). We examined the effect of marker choice on taxonomic identification of prey items and investigated diversity metrics in relation to individual and environmental covariates.
Additionally, we sought to quantitatively compare our results using HTS methods to a Blue Catfish diet study that employed a hybrid approach using visual observations of stomach content coupled with barcoding of single, yet unidentifiable tissues (Schmitt et al., 2018).

| Field collections
In fall 2016 (September-October), we obtained 136 Blue Catfish samples ranging in size from 170 mm to 770 mm TL from two sampling sites located three km apart (hereafter referred to as Chericoke Retreat and Putney's Mill) on the Pamunkey River, VA.
All sampling occurred in tidal-fresh habitats with minimal salinity influence (0-0.1 ppm) at 65 rkm (measured from the York River mouth at the Chesapeake Bay). Each collection transect was approximately 2 km each in length. We used a boat-mounted Smith-Root electrofishing 7.5 GPP unit, 5,000-watt generator at 15 pulses per second. All fish were measured (TL, in mm) and placed into containers with ice to minimize tissue degradation. On the same date, fish were taken to a controlled laboratory setting for GI tract extraction and preservation. Entire GI tracts (esophagus to anus) were carefully removed with sterile procedures, placed in 100% ethanol, and stored in a temperature-controlled facility until further processing.

| Genomic DNA extraction and library preparation
We briefly examined contents from the GI tract and removed them using sterile techniques. Two Blue Catfish contained empty stomachs and digestive tracts and were not included in further analysis. DNA was extracted using the Macherey Nagel Stool Sample Extraction kit. Amplification of barcoding regions was conducted for COI (Leray et al., 2013) and for 18S (Zhan et al., 2014). We employed a species blocking primer 5′-CAAGAATCAGAAAAGGTGTTGGTAA AGA-3′ as outlined in Leray et al. (2013)

| High-throughput DNA sequencing (HTS) analysis
Sequence reads obtained from the MiSeq platform were analyzed using QIIME 2-2019.4 (Bolyen et al., 2019), with 18S and COI sequences split into two separate groups for processing. For each run, primer sequences were trimmed and forward and reverse reads joined and filtered using DADA2 (Callahan et al., 2016). Singletons were removed, and any reads identified in the negative control were subtracted from other samples. Sequences from all six runs were then combined into one table and clustered at 98% using Vsearch (Rognes et al., 2016). Taxonomy was assigned by comparing 18S features against the SILVA 132 database (695,171 sequences, Quast et al., 2013) in QIIME 2, while COI taxonomy was assigned using BLAST + local alignment (Camacho et al., 2009) with a reference database containing 1,769,786 COI sequences downloaded from NCBI.
One feature representing the host species was removed from each sample to eliminate host-introduced bias in the diet analysis.

| Diet diversity and composition
We produced alpha-rarefaction curves for 18S and COI using QIIME 2. For 18S reads, samples were rarefied to 1,355 while COI reads were rarefied to 5,403. Alpha-and beta-diversity analyses were run using a rooted tree in QIIME 2 on rarefied tables. Kruskal-Wallis tests of alpha diversity matrices were performed, and q values (p values with a Benjamini-Hochberg correction, Benjamini & Hochberg, 1995) of less than 0.05 were accepted as statistically significant. For diversity metrics analyzing size, fish were grouped by 100 mm increments.
Beta-diversity comparisons were run with PERMANOVA using a Bray-Curtis and/or a weighted unifrac distance matrix and Adonis as implemented in QIIME 2. Differential abundance was analyzed with ANCOM (Mandal et al., 2015) on nonrarefied tables. Taxonomy bar plots were created using nonrarefied, collapsed taxonomy tables.
We compared the results from our HTS molecular approach to data obtained from a previous hybrid study conducted by Schmitt et al. (2018). This hybrid dataset comprised 617 Blue Catfish diet samples collected in the same general location as our study sites on the Pamunkey River during the months of September and October in 2015. Richness, diversity, and diet composition were compared to HTS results across sampling dates and life stages. Life stage was split by juvenile (up to 300 mm) and adult (over 300 mm). Plots were constructed using the "ggplot2" package in R (Wickham, 2016).
We correlated diet diversity metrics to daily flow measurements (USGS 01673000 gaging station near Hanover, VA). The hybrid study analyzed contents found in stomachs only, which would have been consumed within the past 24 hr. Therefore, we correlated the hybrid dataset to average flow within one day of collections. (Carreon-Martinez et al., 2011). We estimated average flow rate for HTS analysis using the previous three days of each sampling event. We based this decision on a study by Schultz et al. (2015) that demonstrated gut evacuation rate in salmonids was approximately 2 days at 22°C and decreased as temperatures decreased. As the average temperature over our sampling dates was 19.4°C, we felt three days were a more accurate predictor for persistence of diet items within the gut for our samples.
Cumulative prey curves for both genes and approaches were constructed to identify the number of samples needed to characterize Blue Catfish diets. Curves and associated 95% confidence intervals were calculated with EstimateS, version 9.1 (Colwell, 2013), where the cumulative number of unique taxa were plotted against the randomly pooled samples. This process was bootstrapped 1,000 times to generate means and 95% confidence intervals. We used the mean slope (B) of the last five subsamples (linear regression) as an objective criterion for sample size sufficiency, where sample size is considered sufficient when B ≤ 0.05 (Bizzarro et al., 2007;Brown et al., 2012).

| Blue catfish diet analysis
Blue Catfish prey assessment revealed a highly diverse diet that included plant matter, fish, crayfish, turtles, terrestrial insects, aquatic macroinvertebrates, molluscs, algae, and several microscopic phyla ( Figure 2). According to Faith's phylogenetics diversity (PD), fish collected on September 23 or October 20 had more diverse diets than those collected on October 3 or October 12. Faith's PD did not differ between collection sites, size of adult fish, or life stage of fish (juveniles versus adults). Weighted unifrac pairwise tests for both 18S and COI detected differences in diet for the sampling date groupings mentioned above (Table 1). COI sequences also detected distinct diets for collection sites (pseudo-F = 2.64, q = 0.008) as well as for life stage (pseudo-F = 5.30, q = 0.001), with juvenile Blue Catfish consuming noticeably more Asian Clam than their adult counterparts ( Figure 2c).

F I G U R E 2
Blue Catfish diet taxonomic bar plots. Any assigned diet items constituting at least one percent of total diet are shown. Bar plots have been scaled to illustrate the relative proportions of each of these top diet items. Comparison of diet content at phylum (a) and genus (b) levels using either COI, 18S, or the hybrid method illustrates how examined diet content varies depending on methodology. Additionally, COI sequences revealed two ontogenetic shifts in our dataset (c). As Blue Catfish transition from juveniles to adults, we observed a decrease in consumption of Asian Clam. At 500 mm, Blue Catfish began to shift toward piscivory as they decreased plant and mollusc consumption and increased fish and crayfish predation  (Table 3). Collection date in conjunction with collection site could explain an additional 3.9% to 4.5% of the variation observed. COI sequences indicated 3.7% of variation could be attributed to life stage and 6.1% to size of fish.
As collection date consistently appeared to be an important factor   Shad Alosa sapidissima, species of concern under a moratorium in this geographical region due to low population levels (Greene et al., 2009).

| Comparison of HTS and hybrid approach
We compared our taxonomic results to those achieved using a hybrid approach by Schmitt et al. (2018). The hybrid approach identified an average of 1.1 diet items per stomach. In contrast, HTS using DNA extracted from GI tracts resulted in an average of 45.3

F I G U R E 3 (a) Correlation between
Blue Catfish diet richness and flow rates using high-throughput sequencing.
The relationship between flow rates and diet richness was calculated for Faith's phylogenetic diversity, richness as measured by observed operational taxonomic units (OTUs), and Shannon diversity index using both 18S (red) and COI (blue). All analyses demonstrate an inverse relationship where low flow rates are associated with increased diet diversity. (b) Correlation between Blue Catfish diet richness and flow rates using hybrid methodology. The relationship between flow rates and diet richness was calculated using the Shannon diversity index. Despite the lower richness detected with hybrid methods and the lower flow rates observed on sampling dates for this cohort, we observed the same inverse relationship where low flow rates are associated with increased diet diversity as seen with high-throughput sequencing data COI OTUs and 22.7 18S OTUs. Taxonomic analysis indicated that the hybrid method failed to identify microscopic incidental diet items such as ciliates, apicomplexans, and rotifers as well as easily digestible taxa including fungi, sponges, and worms ( Figure 2).
Diets from the hybrid study skew heavily in favor of larger, slower to digest prey such as fish, crabs, and bivalves. Cumulative prey curves were developed for all methods (Figure 4). In order to accurately compare HTS to the hybrid dataset, we removed incidental taxa, including microscopic organisms and parasites, captured by HTS that would not have been the object of direct predation by Blue Catfish and therefore would not have been counted in a traditional morphological study. It was immediately apparent that cumulative prey curves are biased by the taxonomic resolution of the study; that is, tools (like HTS) that more precisely identify prey to species levels will inherently require more samples to achieve Lastly, we examined hybrid diet content in relation to flow rates as performed for HTS results. As seen with HTS data, samples collected on low flow rate days exhibited increased diet diversity when compared to samples collected on higher flow rate days (Figure 3b).
Thus, the negative correlation between Blue Catfish diet diversity and flow rates was observed across multiple years using multiple methods.

| D ISCUSS I ON
Taxonomic analysis of diet contents using HTS revealed greater levels of at-risk species consumption than observed using morphological or hybrid studies. The trophic ecology of invasive Blue Catfish has been investigated in several tidal rivers in Virginia (MacAvoy et al., 2001;Schloesser et al., 2011;Schmitt et al., 2019) and elsewhere in Chesapeake Bay (Aguilar et al., 2017;Iwanowicz et al., 2019) with high spatial and temporal resolution. Despite these efforts, federally endangered Atlantic Sturgeon predation was not detected through morphological identification or DNA barcoding of unidentified diet tissues. Similarly, samples used in the hybrid approach collected from the Pamunkey River did not detect predation of at-risk American shad, Alewife, Blueback Herring, or American Eel Anguilla rostrata.
However, HTS methodology detected considerable consumption of these species. Our study thus supports the hypothesis that traditional morphological and even hybrid approaches may miss important prey, especially small, soft-bodied organisms, eggs, and larvae that digest rapidly (Bromley, 1994;Carreon-Martinez et al., 2011;Legler et al., 2010;Schmitt et al., 2019). Future studies quantifying overall consumption of at-risk species would be beneficial and would require tributary-specific density calculations for both predator and prey species to determine impact levels.

F I G U R E 4
Cumulative prey curves. Cumulative prey curves were constructed for 18S, COI, and the hybrid morphology-based method. As Blue Catfish do not intentionally prey upon incidental taxa and would therefore not have been counted in traditional morphological studies, they were removed from high-throughput sequencing (HTS) data for direct comparison to the hybrid method. Additionally, insects were pooled to family level to match the taxonomic resolution used in the hybrid study. The resulting curve for 18S reached asymptote (B ≤ 0.05) at N = 67, COI at N = 114, and hybrid at N = 204. Both 18S and COI methods identified more unique taxa than the morphology study, yet reached asymptote with fewer sample numbers required, demonstrating that HTS outperforms traditional morphology-based studies Our results suggest that Blue Catfish diet diversity is dependent on daily flow rates, with high flow rates resulting in lower diet diversity. This relationship held true across analysis methods, years, and flow magnitude. We hypothesize that Blue Catfish exhibit higher mobility and prey-seeking behaviors during lower flow in the Pamunkey River. Conversely, they may decrease movement to conserve energy under higher flows, leading to lower probabilities of encountering diverse prey assemblages. Alternatively, higher flow rates could decrease consumption by interfering with Blue Catfish olfactory senses, by providing additional cover for prey through inundation of the aquatic-terrestrial transition zone, or by washing off periphyton film from surfaces such as SAV, thereby decreasing incidental prey consumption.
Strong tidal influence occurs at our sampling sites and other tidal rivers in the Chesapeake Bay. The interplay between daily flow rate and tidal influence coupled with fine-scale feeding ecology provides an opportunity for future research. For example, studies such as weekly feeding chronologies could be analyzed with HTS and correlated to tide schedule over a broad range of flow rates. Additional studies addressing gut evacuation rates in this species and ecosystem would also be of benefit to more precisely correlate flow rates, temperatures, and diet diversity.

Diet diversity analyses indicated multiple ontogenetic shifts as
Blue Catfish size increased. Results indicate that juvenile diets are as diverse as those of adults, yet distinct with Asian Clam abundant in juvenile samples. Similarly, we noted no significant decrease in overall diversity as adult fish grew in length, but beta-diversity analysis did indicate a shift toward piscivory around 500 mm TL. Schmitt et al. (2018) showed an ontogenetic diet shift from omnivory to piscivory between 500 and 900 mm TL depending on river, with Pamunkey River Blue Catfish shifting at 900 mm TL. It is important to note, however, that this shift was determined using N = 4,322 stomachs collected from tidal fresh, oligholaine, and mesohaline waters over the course of four field seasons, whereas the current study was completed at much finer spatiotemporal scales. However, the expanded ability to detect prey items coupled with the finer taxonomic resolution achieved through HTS may allow such ontogenetic shifts to be identified at earlier stages than previously determined. The limited spatiotemporal scope of the current study may also help explain why we observed ontogenetic shifts to piscivory at smaller sizes, as such shifts may be driven by seasonal and(or) spatial trends. Schmitt et al. (2018) defined piscivory as the probability that fish occurred in at least 50% of stomachs examined for a particular size range. However, all size ranges in our study, including our smallest of 100-200 mm TL, met this 50% criteria. These taxa detections may occur from consumption of eggs or larvae and would likely be missed using morphological techniques. As detection of prey life stage is not possible using HTS, these studies may need to use different parameters for defining piscivory, such as percent fish consumed in relation to total diet items.
The lack of decrease in diet diversity as Blue Catfish grew in size is surprising given that, as discussed above, Blue Catfish are generally understood to become less omnivorous and more piscivorous as they grow (Eggleton & Schramm, 2004;Graham, 1999;Schmitt et al., 2019 Despite the larger hybrid diet dataset, we found that HTS provided significantly higher diversity and richness estimates, and ultimately a broader and more holistic view of Blue Catfish diets.
The number of taxa detected increased 20-to 40-fold when using HTS on GI tracts in comparison with hybrid methodology examining stomachs alone. Given the large disparities in diet composition across approaches and the documented abilities of morphologybased methods relative to DNA-based studies, the differences found here are not likely due to annual variation in Blue Catfish diets.
Similar observations have been made for other predator species and ecosystems previously examined (Cowart et al., 2015;Jakubavičiūtė et al., 2017) due to the inability of traditional morphological analyses to identify partially digested items, microscopic organisms, and rapidly digested tissue such as eggs and larvae (Guillerault et al., 2017;Su et al., 2018). Cumulative prey curves demonstrate that HTS requires approximately half the samples needed for morphologybased studies to capture full diet diversity for Blue Catfish.
It is important to note that HTS diversity and richness calculations may be inflated due to the detection of incidental and secondary diet items. For example, SAV contains high biomass of species-rich periphyton (biofilm) on leaf and stem surfaces that contain a diverse assemblage of algae, cyanobacteria, heterotrophic microbes, and detritus (Gordon-Bradley et al., 2014;Hoagland et al., 1982). Thus, SAV provides habitat and food resources for macroinvertebrates, zooplankton, and small fishes. As omnivorous fishes feed on and among these SAV beds, they are also passively consuming incidental organisms within this periphyton assemblage that are likely contributing to the high diet diversity observed with HTS results in comparison with the hybrid approach.
Furthermore, the use of GI tracts in this study, while increasing the ability to detect primary prey consumed over a longer time frame, allows for mixing of primary and secondary prey DNA as items are digested.
However, Jakubavičiūtė et al. (2017) hypothesized that secondary prey contribution to diet diversity would be minimal given their decreased biomass and advanced DNA degradation. Further studies examining the contribution of incidental and secondary prey to predator diets would be of interest as the use of HTS to examine diet diversity increases.
While our study indicates that HTS methods are generally superior to morphological methods given the lower sample sizes required as well as the breadth of results achieved, there may be instances in which traditional analyses are still appropriate. For instance, many diet studies depend on data such as biomass and direct counts (Christensen, 2009;Hyslop, 1980 In addition to preferential amplification of species, the two genes used in this study exhibited differences in taxonomic resolution of diet items. As seen in previous studies (Günther et al., 2018;Holman et al., 2019), the more conserved 18S rRNA gene amplified across a broader range of taxa, leading to greater diversity values when examined at high taxonomic levels in comparison with COI. The broader amplification range of 18S also allowed for higher detection of incidental diet items, which researchers may or may not want to include depending on study objectives. Conversely, the increased sequence variation of COI allowed for a more finite resolution of successfully amplified taxa and ultimately led to a higher number of detected species, despite its narrower phylogenetic scope. The fine-scale taxonomic resolution achieved with COI also detected differences in alpha and beta diversities that were not detected using the more conserved 18S gene.
Our dataset contained many COI sequences that could not be taxonomically assigned to a phylum. These unassigned COI sequences reflect the degree to which the wealth of biodiversity has yet to be genetically catalogued. Mora et al. (2011) hypothesized that only 13.8% of species on earth have been described. Another study estimated that in freshwater systems, approximately 3,000 fish and 100 bivalves remain undescribed (Tedesco et al., 2014). With the advent of metabarcoding, we now have the ability to genetically detect these undescribed creatures. Moving forward, it will be important to curate vouchered specimens with DNA sequences as we seek to increase our knowledge of living organisms and predator-prey interactions.

ACK N OWLED G M ENTS
We

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data available from OSF: https://osf.io/g6wcq. Evans, H., & Bunch, A.
High throughput dual gene sequencing of Pamunkey River (Virginia, USA) tidal-fresh fish diets.