Authors: Mary C. Fisher, Emily W. Grason, Alex Stote, Ryan P. Kelly, Kate Litle, P. Sean McDonald
Categories: Research Article, Biology and Life Sciences, Ecology and Environmental Sciences, Medicine and Health Sciences, Research and Analysis Methods, Biology and life sciences, Research and analysis methods
Source: PLoS ONE
Predation by invasive species can threaten local ecosystems and economies. The European green crab (Carcinus maenas), one of the most widespread marine invasive species, is an effective predator associated with clam and crab population declines outside of its native range. In the U.S. Pacific Northwest, green crab has recently increased in abundance and expanded its distribution, generating concern for estuarine ecosystems and associated aquaculture production. However, regionally-specific information on the trophic impacts of invasive green crab is very limited. We compared the stomach contents of green crabs collected on clam aquaculture beds versus intertidal sloughs in Willapa Bay, Washington, to provide the first in-depth description of European green crab diet at a particularly crucial time for regional management. We first identified putative prey items using DNA metabarcoding of stomach content samples. We compared diet composition across sites using prey presence/absence and an index of species-specific relative abundance. For eight prey species, we also calibrated metabarcoding data to quantitatively compare DNA abundance between prey taxa, and to describe an ‘average’ green crab diet at an intertidal slough versus a clam aquaculture bed. From the stomach contents of 61 green crabs, we identified 54 unique taxa belonging to nine phyla. The stomach contents of crabs collected from clam aquaculture beds were significantly different from the stomach contents of crabs collected at intertidal sloughs. Across all sites, arthropods were the most frequently detected prey, with the native hairy shore crab (Hemigrapsus oregonensis) the single most common prey item. Of the eight species calibrated with a quantitative model, two ecologically-important native species–the sand shrimp (Crangon franciscorum) and the Pacific staghorn sculpin (Leptocottus armatus)–had the highest average DNA abundance when detected in a stomach content sample. In addition to providing timely information on green crab diet, our research demonstrates the novel application of a recently developed model for more quantitative DNA metabarcoding. This represents another step in the ongoing evolution of DNA-based diet analysis towards producing the quantitative data necessary for modeling invasive species impacts.
Impacts of biological invasions can be difficult to measure at the early stages for a variety of reasons. In many cases, managers lack baseline data against which to compare population changes in the local ecological community [1]. Even when such data do exist, per capita impacts may not scale up to statistically detectable community shifts when the overall invasive population size remains relatively small [2]. Additionally, assessment of impacts can “lag” behind their occurrence if impacted species have not yet had enough time to incur detectable demographic shifts (e.g., changes in vital rates; [3, 4]). Yet gauging impacts early in the invasion process can buy managers valuable time in prioritizing limited resources to geographies or habitats before they incur dramatic losses [5].
For predatory invasive species, diet studies offer crucial early insight into the range and magnitude of trophically-mediated community impacts. Diet analysis can also be combined with manipulative experiments and behavioral studies to parameterize population-level impact analysis for later invasion stages. For example, several diet studies of the invasive lionfish (Pterois volitans and P. miles) elucidated the range, electivity, and size-dependency of lionfish diets in the Caribbean ([6–8] and others cited in [9]); when results were combined with behavioral observations and models on diurnal activity patterns [10], prey switching [9] and functional responses [11], researchers were able to predict population-level impacts.
While direct characterization of diet has traditionally relied on visual identification of stomach contents (including [6–8]), DNA metabarcoding is an increasingly popular method to for this purpose. DNA metabarcoding has been shown to identify stomach contents to higher taxonomic resolutions than visual analysis, and can reveal novel predator-prey interactions through improved detection of certain prey groups that are visually indistinguishable after mechanical digestion (e.g., those composed primarily of soft tissue; [12–14]). There have also been recent methodological improvements in the analysis of DNA metabarcoding data more broadly, aimed at the difficulty of translating metabarcoding data–which consists of DNA sequencing read abundances–into accurate quantitative information about the sampled environment. Chief among these challenges is accounting for bias introduced during the required DNA amplification process, which causes sequencing read abundance to be unrepresentative of the true composition of a DNA sample [15, 16]. However, with quantitative models that were recently developed for other applications of DNA metabarcoding, we can now use mock prey “communities” to calibrate sequencing read abundance for such amplification biases, and quantify the contribution of different species’ DNA to a stomach content sample.
We combine DNA metabarcoding of crab stomach contents with quantitative models to describe the diet of invasive European green crabs (Carcinus maenas; hereafter “green crab”) in Willapa Bay, Washington, where there have been recent dramatic increases in crab abundance. Green crabs were first detected in Willapa Bay in 1998 [17], but apparently failed to establish a sustained presence in the estuary during the subsequent two decades [18]. Recent trapping and monitoring efforts, as well as opportunistic observations made by local shellfish growers, detected a dramatic increase in both the relative abundance of green crabs as well as their geographic range within the Bay ca. 2015–2017. These trends have raised concern that shifting conditions are enabling green crabs to transition from a lag phase into population growth, establishment, and spread. These later phases of the invasion process are most commonly associated with impacts to native and commercially important resident species. Regional resource managers anticipate significant impacts of the green crab invasion, with particular concern for the native habitat-building eelgrass (Zostera marina), and shellfish species of high commercial and sociocultural importance. Willapa Bay contributes approximately one quarter of all Washington State shellfish aquaculture production [19], predominantly through the harvest of Manila clam (Ruditapes philippinarum) and Pacific oyster (Magallana gigas; [19]). Shellfish growers have expressed concern that green crab might already be reducing the harvest of natural set Manila clams on which the industry relies.
Though species- and community-level impacts of green crabs have been documented in other parts of their invasive range [20–22], most remain speculative in Washington based on likely or observed habitat overlap between green crab and species impacted elsewhere, or on behavioral experiments that might not translate to population impacts in the field. For instance, McDonald et al. [23] demonstrated that juvenile green crab can both prey on and outcompete size-matched native Dungeness crab (Cancer magister) for food and shelter, but the two species did not co occur at that size during that period. Additionally, while research in California has documented declines in the native hairy shore crab (Hemigrapsus oregonensis) attributable to green crab [20, 21], and the two species frequently co-occur in green crab traps in Washington State [24–26], no baseline data exist to robustly evaluate local green crab impacts to H. oregonensis populations. Since feeding behavior of generalist invasive predators can vary substantially across their global range [27, 28], local, habitat-specific diet analyses are needed to assess the context-dependence of green crab diets and consequent community impacts.
We provide a timely investigation of European green crab diet at clam aquaculture beds and intertidal sloughs in Willapa Bay, Washington, and quantify the consumption of prey species of interest within and across site types. We first identified putative prey taxa from green crab stomach contents using DNA metabarcoding of the mitochondrial cytochrome oxidase I (COI) region. We then constructed mock prey “communities” for a subset of species of particular socioeconomic and ecological concern, and corrected for amplification bias using a novel quantitative model [29] in order to estimate the abundance of each species’ DNA in crab stomach content samples. We used this calibrated abundance data to model the composition of the average diet of a green crab from a clam aquaculture bed versus an intertidal slough. Our work not only provides the first detailed description of in situ green crab diet in Washington State, it also broadens the potential downstream applications of diet DNA (dDNA) to invasive species research and management by incorporating novel quantitative methods from the broader environmental DNA literature.
We collected green crabs at four locations in Willapa Bay, Washington, in May, July, and September of 2021 (Fig 1 and Table 1). Green crab collection was conducted by the Washington Sea Grant Crab Team and the Willapa Grays Harbor Oyster Growers’ Association under Washington Department of Fish and Wildlife Aquatic Invasive Species Permit EGC.001-21.20. Two collection sites were located on tide flats under active cultivation of un-netted Manila clams (R. philippinarum; clam aquaculture bed or “clam bed”), and two were uncultivated intertidal sloughs (“slough” sites) adjacent to emergent salt marsh vegetation. Sites were paired such that one of each type was found on either side of a transition zone within Willapa Bay known as the “fattening line.” Sites to the north of the fattening line (toward the estuary mouth), are characterized by greater influence of ocean water, with lower residence time than those to the south (toward the head of the estuary; [30]). Spring and summer months were sampled under the assumption that temperature-influenced energetic demands increase foraging activity [31, 32], potentially reducing the incidence of empty stomachs and increasing encounter rates with prey of interest.
Fig 1 Map of the four collection sites in Willapa Bay, WA, with the distribution of carapace widths (in millimeters) at each site for all European green crabs trapped (gray), those chosen for sequencing (black), and those chosen for sequencing which had identifiable prey DNA in their stomach contents (red). The map inset shows the location of Willapa Bay on the Washington State coast. The main map background is from the U.S. Geological Survey National Map, compiled using terrainr [33].
Crabs were captured using galvanized steel minnow traps and square Fukui fish traps, baited with Pacific mackerel (Scomber japonicus) and set for an overnight high tide (for detailed trapping protocol, see S1 Text). On trap retrieval, green crabs were sexed and their carapace width measured to the nearest millimeter. They were then immediately placed on dry ice to preserve stomach contents from degradation. Crabs were transported to the University of Washington on dry ice, and stored in a -20°C freezer until dissection. Trap-level bycatch presence and mortality data were also recorded during sampling to support inference of trap-facilitated predation.
Of the 344 total green crabs captured during these efforts, 121 were selected for dissection. We preferentially selected larger, predominantly male crabs for analysis (Fig 1) to maximize the likelihood that crabs would be large enough to handle mature adult Manila clams; green crab prey size has been shown to scale with crab master chela height and carapace width [31, 34]. We removed stomach contents from each crab by inserting a pipette with a filter tip through each crab’s mouth and esophagus, in order to limit the amount of predator DNA and gut bacteria included in the stomach sample for DNA extraction (following [35]). If the stomach was nearly empty, we used gastric lavage (with 400uL 100% EtOH) to help remove particulates. Stomach contents were deposited in a sterilized 1.5mL microcentrifuge tube and stored in 100% EtOH at -20°C until DNA extraction. Stomach fullness was then ranked on a 0 (empty) to 7 (full) scale [36]. All dissection tools were sterilized between individuals with 20% NaClO, 100% EtOH, and flame.
DNA was extracted from stomach contents with the Qiagen DNeasy Blood & Tissue Kit (Qiagen Corp., Valencia, CA, USA), with some modifications to the kit protocol, including an initial drying period in a vacuum centrifuge to evaporate 100% EtOH. After the first PCR step described below, we found that a subset of samples did not yield adequate DNA from the DNEasy Blood & Tissue Kit (i.e., could not be visualized on 1.5% agarose gel or had too low a concentration of amplicon to be quantified in the Quant-iT Assay). These samples were re-extracted using the Qiagen PowerSoil Pro Kit (Qiagen Corp., Valencia, CA, USA).
DNA metabarcoding was completed using a two-step PCR protocol (modified from [37]). The first PCR step used the BF3/BR2 primer pair to amplify a 418bp segment of the mitochondrial Cytochrome Oxidase I (COI) region (within the Folmer fragment; [35, 38, 39]). We chose this primer pair because of its prior application in green crab diet analysis [35], and because preliminary analysis on a subset of our samples showed its ability to amplify DNA from prey species of particular management and conservation interest. The 25uL PCR reaction consisted of 2X Qiagen Multiplex Plus Master Mix, 0.5uM of each primer, and 4uL of DNA extract. After an initial denaturation step (95°C for 5 minutes), we repeated for 35 95°C for 30 seconds; 50°C for 30 seconds; 72°C for 50 seconds). This was followed by a final extension at 68°C for 10 minutes. All DNA extracts were amplified in three separate reactions to produce three technical replicates per crab stomach. PCR product was visualized on a 1.5% agarose gel. Samples with amplicon that produced a visible band on the agarose gel were cleaned to remove primer dimers using Mag-Bind Total Pure NGS beads (Omega Bio-Tek, Norcross, GA, USA), and quantified using the Quant-iT dsDNA High Sensitivity Assay Kit (ThermoFisher Scientific, Waltham, MA, USA). We then diluted all samples with high enough concentrations of cleaned amplicon (>0.89ng/uL) to 10ng DNA in 11.25uL by adding the appropriate amount of UltraPure DNase/RNAse-free water (Fisher Scientific, Pittsburgh, PA, USA). The second PCR added index sequences (6 base pair nucleotide tags; IDT for Illumina DNA/RNA UD Indexes) to the amplicons using 1X HiFi HotStart Ready Mix (KAPA Biosystems, Wilmington, MA, USA) and 11.25uL of diluted (high concentration samples, >0.89ng/uL) or pure (low concentration samples, ≤0.89ng/uL) amplicon. The number of PCR cycles varied between 5, 8, and 12 depending on amplicon concentration. Indexed amplicons were visualized on a 1.5% agarose gel, cleaned, and quantified using the same procedures as the first PCR step. This process yielded at least three laboratory samples per crab of uniquely indexed amplicons. Detailed laboratory protocols are available on GitHub (accessible doi.org/10.5281/zenodo.10967264).
The same two-step PCR process was completed for a PCR positive control (DNA sample from kangaroo, genus Macropus; amplified for n = 9 replicate controls across all separate PCR runs), a PCR negative control (n = 9), and a DNA extraction negative control (n = 7). Positive and negative controls from across DNA extractions / PCR runs were pooled prior to indexing, so that each sequencing run contained one positive control sample and one negative control sample.
We pooled indexed amplicons and controls for 2x300 (paired end) sequencing on an Illumina MiSeq v3. The first sequencing run included indexed amplicons from n = 30 crab (90 total samples) plus controls, and was completed by the University of California, Los Angeles’ Technology Center for Genomics & Bioinformatics. The remaining three sequencing runs consisted of indexed amplicons from n = 56 samples (168 technical replicates) plus controls, and was prepared and sequenced in-house using the MiSeq Reagent Kit v3 600-cycle (Illumina, San Diego, CA, USA) at the Center for Environmental Genomics Laboratory, University of Washington. All pooled libraries sequenced at the University of Washington’s Center for Environmental Genomics were run at a concentration of 7pM with 20% phiX. DNA sequencing data are accessible at doi.org/10.5281/zenodo.10850508.
DNA sequencing data was cleaned and analyzed in R v4.1.3 [40] and Git for Windows v2.35.1 (running bash v4.4.19), using a combination of custom scripts and the programs cutadapt [41], DADA2 [42], and NCBI’s BLAST [43, 44], to translate demultiplexed DNA sequencing data into a filtered dataset of putative prey Amplicon Sequence Variants (ASVs) with corresponding sequence counts (modified from [45]. In short, we trimmed primers from demultiplexed DNA sequences using cutadapt and trimmed sequence length according to base pair read quality with DADA2 in R. We then used DADA2 to filter the trimmed sequencing data with the default parameters (maximum of two expected errors, reads truncated at first occurrence of a quality score of two), and to identify ASVs in the filtered sequences with a parametric error model (built on the sequencing data itself; [42]). We also used DADA2 to identify and remove chimeric ASVs, and then removed ASVs found in the positive or negative controls. We aligned the filtered ASVs to NCBI’s BLAST database using blastn, with a minimum percent identity of 94%, minimum word size of 15, and maximum e-value of 1e-30. We identified the last common ancestor for any ASVs that returned multiple species from blastn.
Additional filtering was completed with custom scripts to remove ASVs aligned to (1) green crab DNA, as we are unable to distinguish between predator DNA and putative prey of the same species; and (2) taxa considered non-target (e.g., bacterial, fungal taxa) or too small to be targeted prey (e.g., diatoms, rotifers, unicellular algae, etc.). We first retained only ASVs that were identified to species, and then added higher taxonomic identifications (to family) back into the data if they were not already represented in the species dataset for a given green crab.
From this processed taxonomic dataset, we modified methods from Duprey et al. [46] to manually categorize each taxon as a (1) native taxon, (2) previously detected non-native taxon, or (3) previously undetected non-native taxon. This classification was performed using multiple databases, the World Register of Marine Species (WoRMS; [47]), the Biodiversity of the Central Coast database ([48]), the National Estuarine and Marine Exotic Species Information System (NEMESIS; [49]), the Encyclopedia of the Puget Sound [50], and the peer-reviewed literature. Non-native taxa were subject to additional filtering, to account for the possible mis-assignment of ASVs representing uncatalogued or poorly represented native taxa to non-native close relatives. For previously detected non-natives: if the percent identity with the BLAST reference sequence was >95%, or if percent identity was <95% but native sister taxa were represented in the BLAST database, the taxon was retained; otherwise, we either retained the genus shared by the matched non-native taxon and any native sister taxa (for species-level identifications), or we removed the taxon from the final dataset (none met either criteria). For previously undetected non-natives: if the percent identity with the BLAST reference sequence was >98% and native sister taxa were represented in the database, the taxon was retained (none met this criteria); otherwise, as above, we either retained the identification to the genus level (n = 3 taxa) or removed the taxon from the final dataset (n = 2 taxa). We provide all ASVs and associated taxonomy removed in each filtering step in S1 Table.
The final filtered dataset consisted of 54 unique taxa associated with 274 total ASVs (accessible doi.org/10.5281/zenodo.10967264). Although the DNA from these taxa could be present in the stomach contents for a number of reasons not limited to targeted predation, which we address in the discussion, we will hereafter refer to these taxa as “prey” that comprise the green crab diet.
We used presence/absence information from the final prey dataset to calculate the α diversity of crab stomach contents per site, and β diversities between sampling sites (using the R package vegan; [51]). When reporting α diversities at the site-level or for the entire dataset, we conservatively removed higher-level taxonomic identifications that may already be represented by species-level identifications (e.g., an identification of Crangon sp. in one crab would not count toward α diversity if any other crab collected from the same site had Crangon franciscorum DNA in its stomach contents). To determine whether diet composition differed significantly among the four sampling sites, we conducted a PERMANOVA and then used a PERMDISP, from the adonis2 function in vegan, to evaluate whether the significance in the PERMANOVA was due to variability in sample dispersion around centroids as opposed to differences in centroid location [51]. We also identified which site pairs had significantly different prey assemblages by conducting post-hoc pairwise comparisons with permutation tests implemented in the R package pairwiseAdonis [52]. We completed all statistical testing first using prey presence/absence and the Jaccard dissimilarity index, and then using the Bray-Curtis distance calculated from the eDNA index, an index of read-count proportions [15, 53]. All code to reproduce these analysis is available on GitHub (accessible doi.org/10.5281/zenodo.10967264).
While the eDNA index allowed us to compare relative sequencing read abundances across sites within an individual taxon, bias introduced during PCR processes means that read abundances cannot be directly compared across different taxa [15, 29]. This limited our interpretation of the relative impact of green crab on different prey items of particular ecological or socioeconomic importance. To overcome these limitations for a subset of prey of interest, we converted proportions of sequenced amplicons (read abundances) to proportions of DNA template (calibrated abundances) using a calibration sequencing run of “mock communities” and a quantitative model developed by Shelton et al. [29]. This model corrects read proportions in DNA metabarcoding data for primer- and species-specific amplification efficiency. Posterior mean estimates from the model represent the proportion of DNA from each prey taxon that was present in the PCR template (i.e., subsampled DNA extract from green crab stomach contents).
The Shelton et al. [29] model requires sequencing data from “mock community” samples that consist of known concentrations of DNA from each species of interest. We constructed mock communities composed of nine the Manila clam (R. philippinarum), soft-shell clam (Mya arenaria), Dungeness crab (C. magister), hairy shore crab (H. oregonensis), sand shrimp (Crangon franciscorum), invasive mud snail (Batillaria attramentaria), eelgrass (Z. marina), shiner perch (Cymatogaster aggregata), Pacific staghorn sculpin (Leptocottus armatus), and prickly sculpin (Cottus asper). These species were chosen because of their frequency in the stomach content DNA metabarcoding data; their frequency as bycatch in green crab traps (representing their potential presence in the green crab prey pool, or subjection to trap-facilitated predation); or because they are of management or conservation interest. Collection of Dungeness crab (C. magister) tissue was completed under Washington Department of Fish and Wildlife Scientific Collection Permit McDonald 21–206; collection of M. arenaria and H. oregonensis was completed with permission from the University of Washington’s Friday Harbor Laboratories. Subsampled, archival tissues from the three vertebrate species were provided by The Burke Museum Ichthyology Collection. We also included predator (green crab) DNA in a subset of mock communities.
DNA was extracted from tissue samples using the Omega Bio-Tek EZ-DNA Mollusc Extraction Kit (invertebrate tissues; Omega Bio-Tek, Norcross, GA, USA) or the Qiagen DNeasy Blood & Tissue Kit (vertebrate tissues; Qiagen Corp., Valencia, CA, USA). DNA extracts were amplified using the same PCR process as for the stomach content samples such that a total of five mock communities of varying compositions (S2 Table) were prepared and sequenced on an Illumina MiSeq v3 at the University of Washington’s Center for Environmental Genomics. Sequencing data for the mock communities was processed using the same bioinformatics pipeline as for stomach content samples (hereafter, “sample data”).
Mock community and sample data were then used to run the Shelton et al. [29] quantitative model in R v4.1.3 [40]. We ran the model using three chains, with 5000 iterations per chain, and the default tree length of 12. The estimated median proportion of DNA attributed to each calibrated taxa, for every crab included in the analysis, is reported in the Supplementary Material (S3 Table). We subset the data by model run according to the total number of PCR cycles used the amplify / index each sample (i.e., replicate); for three crabs which were re-processed at two different total PCR cycles, we used only the data from the replicates run at the lowest total PCR cycles.
The resulting estimated proportions of prey DNA were then used to produce an “average” crab diet for each site type, by fitting a zero-and-one inflated Dirichlet distribution with the R package zoid [54, 55]. The zoid input data consisted of a matrix of estimated read counts for each calibrated prey species for each crab; each set of observations was treated as having been drawn from a multinomial distribution for which the probabilities of the K classes were the mean proportions of DNA per prey taxa estimated by the quantitative model, and the size N was the average read depth across technical replicates. Because of the nature of the data (high variability between individuals, with non-zero observations either very large or very small numbers), the likelihood-based optimizer used by zoid had difficulty finding suitable initial values. We solved this by dividing all read depths in the input matrix by 100, effectively reducing our sample size while facilitating model fitting. This transformation is conservative, in the sense that it does not alter the mean model estimates (i.e., species compositions), but does somewhat inflate credible intervals. We fit separate models for crabs from clam bed sites and crabs from slough sites with the fit_zoid function, using four chains with 10,000 iterations per chain.
We dissected a total of 121 green crabs trapped at our four sampling sites (Fig 1 and Table 1), 114 of which had enough stomach contents to sample; we successfully extracted, amplified, and sequenced DNA from 86 of these stomach content samples, 61 of which contained DNA from putative prey items post-filtering (52 male / 9 female; S1 Fig). Crab stomachs in which we detected putative prey DNA most commonly contained one to three distinct prey taxa (S2 Fig).
We identified a total of 54 unique prey taxa across all crabs sequenced. Approximately 90% of these were identified to the species level, and 9.8% to genus. Arthropoda was the most frequently observed phylum across all sampling sites and months, with 32 crab stomachs (52%) containing arthropod prey (Table 2). Green crabs at the Oysterville clam bed site had the greatest prey α diversity (n = 39 unique prey taxa; Table 1); this was also true when site-level prey α diversity was divided by sample size to account for the fact that α diversity tended to scale with the number of crabs sequenced (S3 Fig). Adjusting for sample size, the Nahcotta (n = 9 taxa; 1.5 average per crab) and Stackpole (n = 16 taxa; 1.1 average per crab) slough sites had the second and third highest prey α diversity.
The single most common prey item was the native hairy shore crab (H. oregonensis), which was present in the stomach contents of crabs from all four sites (Fig 2 and Table 2). After the hairy shore crab, the sand shrimp (C. franciscorum) and a genus of red algae (Neoporphyra sp.) were the most commonly identified taxa across all four sites (Fig 2 and Table 2). At the two clam bed sites, the Pacific staghorn sculpin (L. armatus) and an invasive gammarid amphipod (Grandidierella japonica) were as common as the red algae (Neoporphyra sp.) and the sand shrimp (C. franciscorum; Table 2). We did not detect the Pacific staghorn sculpin (L. armatus) nor the amphipod G. japonica in the stomach contents of crabs from the slough sites. Instead, green crabs trapped at the two slough sites most commonly consumed the hairy shore crab (H. oregonensis), brown (Fucus sp.) and green (Ulva compressa) algae, and a nereid worm (Hediste diadroma; Table 2). Notably, despite trapping over half of the sequenced green crabs on Manila clam aquaculture beds, we only detected Manila clam DNA in one green crab (trapped in September at the Oysterville site).
Fig 2 Frequency of prey occurrence (FO%), summarized (a) for major prey groups (distribution of the frequencies of occurrence for individual taxa within each group, across all sites), and (b) for taxa which were present in four or more crab across all sites. Prey group identification numbers from (a) are annotated at the top of (b). In (b), total frequency of occurrence is broken down by site (bar color).
The traps used to collect green crabs for this study also captured bycatch species, some of which were identified in green crab stomach the Dungeness crab, Pacific staghorn sculpin (L. armatus), hairy shore crab (H. oregonensis), invasive mud snail (B. attramentaria) and saddleback gunnel (Pholis gunnellus). Since bycatch and mortality data were recorded for each individual trap, we can specify when our detection of these species in green crab stomach contents was associated with the co-occurrence of the given green crab / prey species in a trap (S5 Table). All green crabs with Dungeness crab in their stomach contents co-occurred with Dungeness crab in our traps; in contrast, none that consumed the invasive mud snail (B. attramentaria) or the saddleback gunnel (P. gunnellus) co-occurred with those species in the traps. Only one crab of the six which consumed Pacific staghorn sculpin (L. armatus), and four out of the twelve that consumed the hairy shore crab (H. oregonensis), co-occurred with those species in the traps. None of these instances of co-occurrence included observed mortality of the given prey species in a trap.
Patterns in β diversity revealed that green crabs from sites of the same type (clam bed, natural slough) had more shared prey species than green crabs from different site types. β diversity was lowest when comparing the two slough sampling sites, Nahcotta and Stackpole (β = 0.42), and the two clam bed sites, Oysterville and Long Beach (β = 0.47; Table 3). β diversity was higher when comparing diet across site types, with the greatest difference in diet composition between the Nahcotta slough site and the Oysterville clam bed site (0.85; Table 3). In a Principal Coordinates Analysis (PCoA) constructed using prey presence/absence, the two clam bed sites and the two slough sites clustered more closely to each other in coordinate space than to sites of the opposite type (Fig 3A). While there were several strong outliers in the PCoA constructed using the semi-quantitative eDNA index, the general pattern of greater proximity between site centroids within, as opposed to across, site types remained (Fig 3B).
Fig 3 Principal Coordinates Analysis (PCoA) constructed using (a) prey presence/absence and the Jaccard dissimilarity index, and (b) the eDNA index and the Bray-Curtis distance.
Variation in crab diet between clam bed and slough sites was statistically significant. A PERMANOVA, run first with presence/absence data and then with the eDNA index, found significant differences in diet composition between two or more of the sampling sites (p<0.05; Table 4). Post-hoc pairwise comparisons (Table 5 and S6 Table) identified significant differences in diet composition between the crabs trapped at Stackpole (slough site) and those trapped at Long Beach (clam bed site; p < 0.05 for both data sets) and Oysterville (clam bed site; p < 0.05 for presence/absence data set alone). All other pairwise comparisons were non-significant. Our PERMDISPs returned non-significant results (p = 0.054 and p = 0.495 using prey presence/absence and the eDNA index, respectively; S7 Table), verifying that the significant PERMANOVAs were not due to differences in the dispersion of samples around the centroid of the PCoA (Fig 3).
To quantitatively compare how different prey of interest contributed to green crab diets, we calibrated sequencing read abundances to account for PCR amplification bias with our mock community data and the Shelton et al. [29] quantitative model. This mechanistic model calibrates sample data relative to a single reference species [29]; a higher amplification efficiency means that a given species may be over-represented in the sequencing read data compared to the reference. With Manila clam as the reference species, green crab had the highest relative amplification efficiency for our BF3/BR2 primer set, followed by the sand shrimp C. franciscorum and Dungeness crab (C. magister; Table 6). The native eelgrass Z. marina was the only species with a lower amplification efficiency than Manila clam (Table 6).
Of the 61 green crabs with prey taxa in their stomach contents, 34 contained DNA from one or more of the calibrated prey species (i.e., those for which we had amplification efficiencies). We did not detect the prickly sculpin (C. asper) in green crab stomach contents, and although the family Zosteraceae was present in stomach contents, it would be inappropriate to apply our calibration derived from Z. marina to an identification that may represent a mixture of different eelgrass and surf grass species (e.g., Z. japonica or Phyllospadix torreyi, both of which are present in Willapa Bay). On average, the eight calibrated prey species present in green crab stomach contents made up less than 50% of the unique prey taxa found in each crab stomach, but there was wide variation in this statistic among crabs; for some crabs collected at the two clam bed sites, the only putative prey items detected in their stomach contents were those which we had included in the mock communities for calibration (Fig 4C). Among crabs collected on clam beds, the calibrated prey represented, on average, almost 100% of all sequencing reads from putative prey taxa, whereas for crabs collected at the slough sites, 50% of reads on average came from non-calibrated prey taxa (Fig 4C).
Fig 4 The diet of an “average” green crab which contains the eight calibrated prey species, from each site type (a). Proportions for each prey species are provided in white text in the appropriate block, and represent the mean mixture proportions estimated from the fitted zero-inflated Dirichlet regression model. To put into perspective what subset of all sequencing data is represented by the eight calibrated prey species, we also provide (b) the proportion of sequencing reads remaining after each filtering step, out of the total aligned to NCBI’s BLAST database (“total reads”), and (c) the proportions of all unique prey taxa, and of reads belonging to all prey taxa, captured by the eight calibrated prey species. 95% credible intervals associated with the mean values in (a) can be found in S6 Fig.
As in the full dataset, the hairy shore crab (H. oregonensis) was the most frequently occurring prey item of those calibrated with the Shelton et al. [29] model (n = 11 / 34 green crabs), followed by the sand shrimp (C. franciscorum; n = 9 / 34 green crabs). When it was detected in a green crab stomach, the sand shrimp also contributed, on average, the highest proportion of DNA out of all calibrated prey species, followed by the Pacific staghorn sculpin (L. armatus) and the Manila clam (S4 Fig and S8 Table). However, there was high variability in calibrated DNA proportions across crabs, and in sample sizes across prey species (S4 Fig).
The sand shrimp (C. franciscorum) and the Pacific staghorn sculpin (L. armatus) were similarly ranked in the calibrated DNA data and the uncalibrated sequencing read data; these species had the two highest proportions of sequencing reads / calibrated DNA, on average, when they were detected in crab stomach contents (S8 Table). In contrast, the Manila clam increased in rank between the uncalibrated (rank = 8/8) and the calibrated (rank = 3/8) data set, while the rank abundance of the soft-shell clam (M. arenaria) decreased from the uncalibrated (4/8) to the calibrated (8/8) data set due to its high relative amplification efficiency (S4 Fig and S8 Table). The difference between a given taxon’s proportion of sequencing reads versus calibrated DNA in each stomach content sample tended to be below 25%, with a few exceptions for which the difference was > 50% (associated with particularly high or low relative amplification efficiencies; S5 Fig).
Modeled average diet composition at the clam bed sites showed nearly even contributions to the diet from each of the eight calibrated prey species (Fig 4A and S3 Table), with the sand shrimp (C. franciscorum) contributing a slightly higher mean proportion of DNA than the other species (mean = 0.14, 95% credible interval (0.12,0.18)). In contrast to the clam bed sites, only three of the eight calibrated prey were identified in the stomach contents of green crabs collected at the two slough sites. The majority of the DNA in an “average” crab stomach from the slough sites was estimated to come from the native hairy shore crab (H. oregonensis; mean = 0.42, 95% credible interval (0.31,0.54)), followed by the sand shrimp (C. franciscorum; mean = 0.32, 95% credible interval (0.32,0.24)), and the soft-shell clam (M. arenaria; mean = 0.25, 95% credible interval (0.16,0.35); Fig 4A and S3 Table).
Our research provides the first detailed diet information for the invasive European green crab in the northeastern Pacific. These data are particularly valuable in Washington State, where an “exponential” increase in green crab populations prompted an emergency proclamation by the state governor (WA Governor’s Emergency Proclamation No. 22–02 Jan 19 2022). We found that green crab in Willapa Bay Washington have a generalist, omnivorous diet that can be highly variable across individuals and locations. Even the most common taxon identified in the stomach contents, the native hairy shore crab (H. oregonensis), had a frequency of occurrence of only ~ 20%. The role of green crab as a generalist predator in Willapa Bay agrees with existing research from other locations where the species is invasive [35, 56], and within its native range [57, 58]. However, our use of DNA metabarcoding allowed for a greater number of prey taxa to be identified from a smaller number of crabs [57] and at a higher taxonomic resolution [56] than prior research using visual diet analysis.
We also distinguished statistically significant differences in stomach contents between sampling sites; green crabs trapped in the natural slough at the Stackpole site had significantly different diet compositions than those trapped on Manila clam aquaculture beds at Long Beach and Oysterville. Although we did not concurrently examine ecological communities at sampling sites, we suspect this is reflective of different community compositions on Manila clam beds and uncultured slough sites, rather than indicative of differing prey selectivity among green crabs at these sites. Previous work on green crab diet on the Patagonian coast found strong similarities between their most frequently occurring prey items and surveys of the local ecological community during separate benthic monitoring [35].
The only other study of green crab diet using DNA metabarcoding (diet DNA or dDNA) was conducted on the Patagonian Atlantic coast [35]. As in the present study, Cordone et al. [35] reported Arthropoda as the most frequently identified phylum in green crab stomach contents, and noted the diversity of algal taxa consumed, particularly the class Phaeophyceae, which had not been documented in prior visual diet analyses. Cordone et al. [35] suggested that algae in the stomach contents most likely represents incidental ingestion, based on their concurrent findings from stable isotope analysis that the energy pathway for macroalgae was less important than that for phytoplankton. However, this does not necessarily exclude the potential for green crab to have negative impacts on macroalgal abundance and diversity.
In contrast to Cordone et al. [35], bony fishes had relatively high frequencies of occurrence in our dataset, particularly the Pacific staghorn sculpin (L. armatus). Other studies have found juvenile fish in green crab stomach contents [57, 59], with increasing crab sizes associated with more frequent fish consumption [59]. Our sampling process, in which we selectively analyzed larger crabs, may therefore explain the greater relative importance of fish species in our dataset than from Cordone et al. [35] (14.5% FO, compared to 5.71% FO). Yet where Taylor [59] estimated that green crab predation has a limited impact on winter flounder (Pseudopleuronectes americanus) populations in Connecticut, there may be greater potential for population-level impacts in Willapa Bay, given the frequency of L. armatus consumption. Moreover, the authors have personally observed Pacific staghorn sculpin aggregating near actively foraging crab, and believe it is possible sculpin may be opportunistically captured and consumed while crabs are foraging for other prey.
The present work was initiated because of concern for Manila clam aquaculture in Willapa Bay, where unexplained declines in production have occurred. Previous work in California has shown green crab to be effective predators on cultured Manila clams [60], and recent increases in the abundance of green crab in coastal Washington waters have been hypothesized to contribute to the observed clam production declines. However, the present study does not provide substantial support for this speculation; only one sampled green crab showed evidence of consuming Manila clam. Notably, this observation occurred in an actively cultivated clam aquaculture bed (Oysterville). Both clam beds sampled were under active culture with a wide size range of clams present, and no netting or bags were used to protect clams from predatory crabs (as is often done in other aquaculture settings).
It is possible that green crabs could be exerting higher predation rates on clam beds outside of our summer sampling period; though consumption and activity rates of green crabs drop in colder temperatures [31, 32], their migration to lower elevations often means that larger relative abundances are observed on Willapa Bay Manila clam beds in winter months. Additionally, size class matches between predator and prey remain an open question for this system. We selected large crabs for diet analysis, but small Manila clams may only be profitable for smaller size classes of green crabs; thus impacts by green crabs may be stage and context dependent, and most severe in habitats where smaller green crabs are found on commercial clam beds. As a result, while Manila clams were an extremely rare diet item in the green crabs we analyzed with the BF3/BR2 primer set, more extensive sampling throughout the year and across smaller size classes may yield different results. Manila clam did have a lower relative amplification efficiency of all but one of the species calibrated in the Shelton et al. [29] model, particularly compared to green crab itself; yet we had no missed detections and adequate read depths of Manila clam in mock communities, which included the species’ DNA at low relative proportions to green crab and other species (e.g., Community D; S2 Table). Given the strong signal of Manila clam in mock communities, we have confidence that if Manila clam had been recently consumed by the green crab we analyzed, we would have detected it in the stomach contents.
In contrast to Manila clam, we detected soft-shell clams (M. arenaria) in the diet of multiple green crabs across several sites. Green crabs have been implicated in the collapse of a soft-shell clam (M. arenaria) fishery in Maine in the mid-twentieth century, when clam abundance declined by 50% in a 4-year period [22]. Green crabs have been shown to preferentially consume soft-shell clams in some studies [61, 62] and substantial effort has been allocated toward strategies to reduce predation in the northeastern Atlantic [63, 64]. Soft-shell clams are non-indigenous to the northeastern Pacific, but have established large populations across Willapa Bay since their initial introduction prior to 1880; soft-shell clams may therefore represent an important food resource for green crab in this estuary.
Although we did detect three cases of green crab consumption of Dungeness crab (C. magister; at Long Beach and Oysterville, in May), Dungeness crabs were recorded as bycatch in the same traps corresponding to these three green crabs. It is therefore possible that consumption of Dungeness crab (C. magister) was facilitated by trapping, and predation would not otherwise have been observed. Notably, however, we did not find physical evidence of predation occurring within the crab traps (i.e., carapace fragments or other debris). We also limited soak times and retrieved traps promptly to reduce the likelihood of consumption occurring within the traps. Nevertheless, hand capture could reduce the likelihood of diet observations being confounded by predator-prey interactions attributable to sampling technique.
Other arthropods observed in stomach contents from this study, while not a focus of commercial aquaculture and fisheries, play an important ecosystem role in Pacific Northwest estuaries. In the present study, a genus of sand shrimps (Crangon spp.) occurred in more than 20% of all samples–primarily as C. franciscorum, when identified to the species level. Our quantitative modeling results further show that when the sand shrimp (C. franciscorum) was consumed by a green crab, the species comprised, on average, a high proportion of DNA in the stomach content sample (compared to other calibrated prey). This particular species (C. franciscorum) is the predominant sand shrimp in Willapa Bay and an important predator, feeding on amphipods, bivalves, and small or juvenile fish [65]. The species is also an important food resource for Dungeness crab (C. magister [36]), sturgeon, and other fishes. Their sister taxon, C. crangon, is known to be an important diet component in the green crab’s native European range [58] in all seasons. We also found that a diversity of amphipods occurred in 29% of the sampled green crabs across all sites. As a taxonomic group, amphipods play a critical role in nutrient cycling and other ecosystem functions and are important prey to a wide variety of invertebrates, fishes and birds. Green crab is known to frequently consume amphipods in both the native [66] and invasive range [56, 67], and significantly reduce their densities in enclosure experiments [68].
Prior work has also demonstrated the negative impacts of green crab on native shore crabs, one species of which (H. oregonensis) was the most frequently occurring taxon in our dataset of green crab diets. In the northeastern Pacific, impacts to native shore crab have been documented in central California [20, 21]. In particular, native shore crab numbers declined and distribution shifted into the high intertidal zone as green crab became more abundant. Moreover, the average size of the native shore crab declined at the same time [21]. Given this evidence of green crab predation on H. oregonensis. in other locations on the U.S. West Coast, as well as the high frequency of H. oregonensis in Washington State green crab traps [24–26], it seems likely that similar population-level impacts could also occur in Washington State estuaries.
One difficulty in interpreting diet data is that we cannot determine whether the presence of a species in the stomach contents is the result of direct or intentional predation. For instance, the detection of some taxa might occur because they were attached to (epibiont), inside of (themselves prey of green crab prey), or otherwise associated with items green crabs targeted as prey. For example, Ansell et al. [57] dissected a complete juvenile fish in a green crab stomach, and found that the fish had previously consumed copepods, ostracods, and an amphipod. We detected all three of these taxonomic groups in our green crab stomach contents. However, we developed our post-BLAST filtering process to select for taxa which were most likely in the size range of true prey for the crabs we sampled, and previous work with terrestrial arthropods suggests that using amplicons greater than 200bp in length reduces the likelihood of detecting DNA from scavenged prey [69], which may have levels of DNA degradation similar to items representing secondary predation. Incidental ingestion might also account for some detections in our study, our identification of Zosteraceae (a family of surfgrasses and eelgrasses) DNA in particular. While green crabs are known to eat eelgrass, both as seeds and blades in experimental [70] and field enclosures [71], respectively, several authors have suggested that consumption may occur incidentally during other foraging activities [72–74]. As mentioned above, algae may also be accidentally consumed during other foraging activities; generally this category represents a small component of the diet relative to known preferred prey like molluscs [35, 56, 58].
Although beyond the purpose of our study, which focused on green crab predation on other organisms, one key aspect of crab diet that cannot be captured using DNA metabarcoding is cannibalism. Cannibalism is common in decapod crustaceans, and has been shown to cause high density-dependent mortality in green crab [75, 76]. However, DNA metabarcoding cannot distinguish between individuals of the same species, and so all green crab DNA in our data was assumed to originate from the predator individual. Future research with the goal of understanding cannibalism in Pacific Northwest green crab populations could employ other molecular techniques to do so, such as genotyping of microsatellite loci, which has revealed cannibalism in invasive lionfish populations [77].
Calibrating DNA metabarcoding data to account for species-specific primer amplification bias provided key insights into our own data set, and contributes important information for future studies of green crab predation. The rank abundance of prey species’ contributions to green crab stomach contents changed after we calibrated sequencing read abundances; both the hairy shore crab (H. oregonensis) and the Manila clam increased in rank abundance in the calibrated DNA data compared to the sequencing read data, whereas the soft-shell clam M. arenaria decreased in rank abundance (S8 Table). Our calibrated data set revealed that out of the eight calibrated species, the sand shrimp (C. franciscorum), the Pacific staghorn sculpin (L. armatus), and the Manila clam contributed the highest average proportion of DNA when present in the stomach contents. By applying the calibrated data to model an average green crab stomach sample, our research also emphasized differences in relative abundance of eight key prey items between Manila clam beds versus natural slough sites; at slough sites, the native hairy shore crab (H. oregonensis) had a higher average abundance in the stomach content DNA than the soft-shell clam (M. arenaria), whereas at the clam bed sites, we did not observe clear differences in average DNA abundance between the eight calibrated prey species.
In generating the calibrated data set, we also revealed that green crab itself had the highest PCR amplification efficiency of all the species included in the mock communities; this was not entirely unexpected, as the BF3/BR2 primer pair was initially developed to study benthic invertebrates [78], but it does suggest that exploring alternative primer sets may help future work minimize the loss of sequencing reads to predator DNA. With mock communities and the Shelton et al. [29] model, future studies can use different primers and still be comparable to this work. The calibration of additional prey taxa that were particularly prevalent in the diet of crabs collected for this study–such as the red algae genus Neoporphyra, and the invasive amphipod G. *japonica–*will be important for assessing the broader ecological impact of green crabs in the Pacific Northwest, and will help put into perspective the relative contribution of the prey calibrated for this study in DNA extracts from green crab stomachs.
A major target for applying diet DNA to invasive species research and management is population-level impact analysis, and the potential for DNA metabarcoding to provide quantitative consumption data for modeling methods like bioenergetics. As the first published application of the Shelton et al. [29] model to stomach content data, our research brings the study of invasive European green crab one step closer to this goal; however, the proportion of observed DNA contributed by each prey species (the “calibrated DNA abundance” output from the Shelton et al. [29] model) is not yet translatable to the proportion of consumed biomass contributed by each prey species (input for bioenergetic modeling). Without accounting for changes in detectability as prey DNA degrades during digestion, quantitative interpretations of dDNA may overemphasize predation impact on prey with lower DNA digestion rates, and, conversely, underestimate predation impact on prey with higher DNA digestion rates [69]. A meta-analysis of spiders found that accounting for differential prey detectability could change dietary proportions of a given spider prey item by up to ~10% [79]. Experimental feeding trials that explore how the detectability of prey DNA changes with time since ingestion, often quantified using the genetic or detectability half-life [69, 80], will be key to additional applications of our research and the broader integration of the Shelton et al. [29] model into invasive species impact analysis. Currently, there are no published feeding experiments to quantify prey detectability half-lives for green crab. Since detectability half-life varies across prey species, feeding experiments to determine prey half-lives can be particularly prohibitive for generalist opportunists like green crab; however, a recent meta-analysis for terrestrial arthropods suggests that there may be consistent and predictable effects of prey traits and environmental conditions on detectability half-lives that make it possible to forego exhaustive feeding trials leveraging meta-analyses and quantitative models [79].
We would like to thank Heckes’ Clam Company, Taylor Shellfish, and Kim Patten for providing property access for sample collection in Willapa Bay, and David Beugli and the Willapa Grays Harbor Oyster Growers Association for assisting with green crab collection. Thank you also to Megan Dethier at Friday Harbor Laboratories, Katherine Maslenikov at the University of Washington’s Burke Museum Fish Collection, Roger Fuller at the Padilla Bay National Estuarine Research Reserve, and Bryan Briones Ortiz at the University of Washington, for their assistance obtaining individuals, tissue samples, and DNA required for the mock communities. Bettina Thalinger, Georgina Cordone, Eily Allen, Erin D’Agnese, Megan Shaffer, and Maya Garber-Yonts shared invaluable lab and bioinformatic knowledge and assistance, and Eric Ward provided help troubleshooting the R package zoid.
Raw sequencing data are publicly available and archived on Zenodo (DOI:10.5281/zenodo.10850507). Code, laboratory protocols, and processed data sets are publicly available on Github, in a release of a repository www.github.com/mfisher5/Green-Crab-dDNA.
This project was supported by the Washington Department of Fish and Wildlife (No.21-18297) as part of the Willapa Bay Oyster Reserves Grant Program (RFP No. 21-0009): https://wdfw.wa.gov/. Mary Fisher was also supported by the NSF’s Graduate Research Fellowship Program (Grant DGE-1762114) while conducting this https://www.nsfgrfp.org/ The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Raw sequencing data are publicly available and archived on Zenodo (DOI:10.5281/zenodo.10850507). Code, laboratory protocols, and processed data sets are publicly available on Github, in a release of a repository www.github.com/mfisher5/Green-Crab-dDNA.