During animal migration, ephemeral communities of taxa at all trophic levels co-occur over space and time. The interactions between predators and prey along migration corridors are ecologically and evolutionarily significant. However, these interactions remain understudied in terrestrial systems and warrant further investigations using novel approaches. We investigated the predator-prey interactions between a migrating avivorous predator and ephemeral avian prey community in the fall migration season. We tested for associations between avian traits and prey selection and hypothesized that prey traits (i.e., relative size, flocking behavior, habitat, migration tendency, and availability) would influence prey selection by a sexually dimorphic raptor on migration. To document prey consumption, we sampled trace prey DNA from beaks and talons of migrating sharp-shinned hawks Accipiter striatus (n=588). We determined prey availability in the ephemeral avian community by extracting weekly abundance indices from eBird Status and Trends data. We used discrete choice models to assess prey selection and visualized frequency of prey in diet and availability on the landscape over the fall migration season. Using eDNA metabarcoding, we detected prey species on 94.1% of the hawks sampled (n=525/588) comprising 1396 prey species detections from 65 prey species. Prey frequency in diet and eBird relative abundance of prey species were correlated over the migration season for top selected prey species, suggesting prey availability is an important component of raptor-songbird interactions during fall. Prey size, flocking behavior, and non-breeding habitat association were prey traits that significantly influenced predator choice. We found differences between female and male hawk prey selection, suggesting that sexual size dimorphism has led to distinct foraging strategies on migration. This research integrated field data collected by a volunteer-powered raptor migration monitoring station and public-generated data from eBird to reveal elusive predator-prey dynamics occurring in an ephemeral raptor-songbird community during fall migration. Understanding dynamic raptor-songbird interactions along migration routes remains a relatively unexplored frontier in animal ecology and is necessary for conservation and management efforts of migratory and resident communities.

Predator Diet Sampling

The Marin Headlands hosts the most significant raptor migration flight along the Pacific Coast of North America where migrating raptors funnel and converge before gaining altitude to cross the San Francisco Bay (Goodrich & Smith, 2008). We collected diet data from raptors banded at a raptor migration monitoring station located in the Marin Headlands operated by the Golden Gate Raptor Observatory/Golden Gate National Parks Conservancy in cooperation with the US National Park Service (37.8262°N, 122.4997°W; Fig. 1). Sharp-shinned hawks were lured into the sampling site with rock doves Columba livia, European starlings Sturnus vulgaris, and house sparrows Passer domesticus, then captured in dho-ghazzas, mist-nets, or bownets, and released after banding and processing (Hull & Bloom, 2001).

We collected eDNA samples from the exterior surfaces of beaks and talons of sharp-shinned hawks (n= 588; Fig. 2) during fall migration in 2015 (n=282) and 2016 (n=276) from September through November. To sample prey DNA, we moistened a sterile histobrush (#25‐2188 Puritan Medical Products Company) in 0.7ml ultrapure water contained in 1.5ml pop-top centrifuge tube. Next, we gently swabbed the entire exterior surface of the upper and lower mandible, targeting any visible prey blood or tissue if present. We avoided contact with interior mouthparts and saliva to minimize collection of predator DNA. We then swabbed the entire surface of each talon, targeting visible prey blood, tissue, or feathers if present. Toe pads or scales were additionally swabbed if visible prey remains were present. After swabbing the raptor, we cut off the nylon brush tip into a 1.5ml screw‐top centrifuge tube containing 0.7ml Longmire's lysis buffer (100mM Tris pH8.0, 100mM EDTA, 10mM NaCl, 0.5% SDS, 0.2% sodium azide) or Queen’s lysis buffer (100mM Tris, 100mM NaCl, 100mM sodium EDTA and 10% n-Lauroylsarcosine, pH 8.0; Suetin et al., 1991) and stored at −20°C. Sampling occurred even if beaks and talons appeared clean (Bourbour et al., 2019, 2021).

We conducted all aspects of this research in accordance with Institutional Animal Care and Use Committee (IACUC; permit #: CA_GOGA_Ely_Raptors_2020.A3), California Department of Fish & Wildlife (California State Permit #: SCP 13739), National Parks Service (study #: GOGA-00004; permit #: GOGA-2022-SCI-0019), and United States Geological Service guidelines (federal bird banding permit #: 21827).

Prey DNA Extraction, Amplification, and Sequencing

We used QIAmp DNA Mini Kit (QIAGEN Inc.) to extract prey DNA from swab tips. We conducted lab work in house at the UC Davis Genomic Variation Laboratory genetics lab that does not process songbird DNA to minimize risk of contamination. We targeted a 464-base pair (bp) amplicon region of the cytochrome c oxidase subunit I (COI) gene using primers COI-fsdF and COI-fsdR which have been demonstrated to have a high-resolution power for identifying species across Passeriformes and Columbiformes (González‐Varo et al., 2014; Bourbour et al., 2021). We modified the primers to have an overhang sequence that would anneal to indexed Illumina adapters (Illumina, 2013; Supplementary Materials Table 3). We tested primers using avian tissue samples from the Museum of Wildlife & Fish Biology at UC Davis. We used orange-crowned warbler Vermivora celata and Swainson’s thrush Catharus ustulatus DNA as positive controls alongside negative controls (PCR-grade water) during library preparation to confirm detection of probable prey species and assess potential contamination or misassignment of reads. We followed the two-step PCR amplification protocol outlined in Illumina (2013), see Supplementary Materials 1 for detailed PCR protocols. After library preparation, we quantified DNA using Quant-iT PicoGreen dsDNA Reagent (Thermo Fisher Scientific) with an FLx800 Fluorescence Reader (BioTek Instruments). We then normalized each sample post-Index PCR to 5ng/µl individually and pooled two libraries at a concentration of 5nM in 15µl of EB buffer (10mM TRIS; pH=8.0-8.4). We sequenced the libraries on two lanes using Illumina’s MiSeq PE300 (v3) platform.

Reference Library and Bioinformatics

Developing custom reference sequence libraries can create more robust results by reducing overrepresented species and representing cryptic diversity (Elbrecht & Leese, 2017; Macheriotou et al., 2019). To start, we compiled a custom reference library of bird species (n=205) that range in the Pacific Flyway according to species account range maps (Billerman et al., 2022; see Supplementary Materials Tables 1,5). We used the R package PrimerMiner (Elbrecht & Leese, 2017) to download all publicly available COI barcode and mitochondrial genome sequences from NCBI and BOLD databases to take full advantage of full and partial sequences available for each species. Sequences were clustered into operational taxonomic units using a 3% sequence similarity to reduce overrepresentation of sequences while capturing sequence variants representing cryptic diversity (Elbrecht & Leese, 2017). Out of 205 avian species, we were able to compile 199 species’ barcode sequences for our reference library. We then manually reformatted the datafiles to be compatible with the reference database format used by the R package DADA2 (Callahan et al., 2016; Supplementary Materials 2).

We filtered out low quality scores (<30) and reads below 250bp using the program Cutadapt (Martin, 2011) and used DADA2 to filter out samples with >2 erroneous base calls, remove chimeras, and merge forward and reverse reads. We matched all barcode sequences to our custom reference library with >99% bootstrap support using the ‘assignTaxonomy’ command in DADA2. Based on assessment of positive and negative controls, we removed prey species detections with <100 total assigned reads and we removed prey species if they represented 0.5% or less of the total number of reads in each individual sample.

Prey availability

To obtain an index of weekly prey availability, we extracted abundances of prey species from the eBird Status and Trends data from August–November using the R package ebirdst (Fink et al., 2020; Strimas-Mackey et al., 2021; eBird Status and Trend 2005 to 2020). The persistence of eDNA in the environment is variable and system dependent (e.g., Andruszkiewicz et al., 2017; Barnes et al., 2014; Strickler, 2015) and more studies are needed to determine the rate of degradation of eDNA in terrestrial ecosystems (Beng & Corlett, 2020). However, because prey DNA on beaks and talons represents diet from previous meals and may be detectable for multiple days on a migrating raptor (Bourbour et al., 2019; Bourbour et al., 2021; Valentin et al., 2021), prey abundance data was extracted from several counties north of the sampling site where sharp-shinned hawks likely occurred prior to their capture at the monitoring station as they follow the coastal mountains (Marin, Napa, Sonoma, Lake, and Mendocino; Fig. 1; Goodrich & Smith, 2008). For this process, we confined eBird data extraction in 2.96km2 spatial cells (which is the native resolution from eBird Status and Trends products; Fink et al., 2020) that also contained sharp-shinned hawk occurrence data according to eBird checklists in the defined study region. Next, we summed weekly relative abundances of prey species across those spatial cells to represent an index of prey availability in the study area.

Statistical analyses

We performed all statistical analyses using R version 4.1.0 (R Core Team, 2021) in RStudio version 2022.2.3 (RStudio Team, 2022). Because the eBird Status and Trends data uses data from 2005-2020 (Fink et al., 2020), we combined both diet sampling years together. We excluded European starling and house sparrow detections from statistical analyses because we cannot confidently rule out contamination at the sampling site as the cause of their detection (rock dove DNA was not detected; Bourbour et al., 2021). We calculated rarefaction and Coverage-based R/E (Type 3) curves with a 95% confidence interval using the R package iNext (Chao et al., 2014) to assess sampling effort. Because sharp-shinned hawk females are 30-40% larger than males (Bildstein et al., 2020), we first tested for differential prey size selection. We used a linear mixed-effects model with average prey mass (extracted from Tobias et al., 2022) as the dependent variable, sex as a categorical explanatory variable, and sample ID (individual hawk) as a random effect using the R packages lme4 (Bates et al., 2015), lmerTest (Kuznetsova et al., 2017), and afex (Singmann et al., 2022), and visualized the model using the R packages effects (Fox & Weisberg, 2019).

Prey traits

We categorized each species detected according to relative size, flocking behavior, and migratory tendency (Supplementary Table 1). For prey size, we used average mass (extracted from Tobias et al., 2022) to classify each species as small (<30g), medium (between 30g and 60 g), or large (>60g). For non-breeding flocking behavior, we used Billerman et al., (2022) to classify species as gregarious if they are described as highly social with conspecifics throughout fall migration, aggregate if described as joining small flocks during the non-breeding season, or solitary if the species is described as solitary and/or seldomly joining mixed-species flocks during migration. For migratory tendency, we used Billerman et al., (2022) to classify species as resident or migratory (i.e., which included both complete and partial migrants) in the Pacific Flyway.

Prey selection

We fitted separate models for male and female hawks given distinct body size differences and statistically significant differential prey size selection. For the models for each sex, we used a set of multinomial logistic regressions (discrete choice models) using the R package mlogit (Croissant, 2013). For each model, we computed the variance-covariance matrix of the parameters to account for repeated measures using the R package sandwich (Zeileis, 2006). These models were used to predict the probability that a species was detected on the migrating predator as a function of species availability and traits. We used prey detection as the dependent variable and included the following as explanatory variables: index of prey availability (from eBird Status and Trends Data), non-breeding flocking behavior, non-breeding habitat association, and migratory tendency. To visualize the diet composition and co-occurrence of selected prey species, we plotted the weekly relative abundances of prey species selected by both males and females over time along with the weekly proportion detected in the diet. We used the scale function in R to normalize (z-score) both eBird relative abundances and proportions in diet. Using the z-scores, we then calculated the Pearson’s R for each of the top seven prey species to measure the correlation between the weekly proportions in diet and eBird prey availability indices (Freedman et al., 2007).