In total, we collected 20 nests from six different bird species (Table 1). Of those, 16 were in nest boxes and four were found over ledges in open areas (American Robin: Turdus migratorius Linnaeus, 1766 and Eastern Phoebe: Sayornis phoebe Latham, 1790). Nests were collected between June 12 and June 28, 2019, as well as on October 5, 2019 at the University of Guelph Arboretum (Guelph, Ontario, Canada, 43.541299, -80.214540) after the fledglings had left or at the end of the breeding period, respectively. Additionally, three abandoned nests were collected in Guelph in July 2018, (Table 1). Nests were individually stored in plastic bags and transported back to the lab, where nest not going into emergence traps were immediately frozen at -20 °C.

Nests collected during June 2019, were transferred into an emergence trap combined with a Berlese-Tullgren funnel (Fig. 1) and left there for a period of 2–3 weeks to allow emergence of arthropods. Contrary to most Berlese-Tullgren funnel set ups, no extra light source was used in order to slow the arthropods’ downward movement and to allow for more time for emergence. Arthropods were collected both in the Berlese-Tullgren funnel and in the collecting bottle of the emergence trap. They were transferred into 95% ethanol for preservation (Fig. 2). Both samples were combined per nest for subsequent DNA barcoding. After the emergence period, each nest was transferred into a fresh plastic bag and frozen at -20 °C until further processing in the laboratory.

Figure 1.

Emergence trap with Berlese-Tullgren funnel (left picture) and nest separated into components >5 mm, >2 mm and <2 mm via consecutive sieving. Individually removed arthropod remains are not displayed.

Figure 2.

Workflow describing the processing of the bird nests for emergence trapping and molecular analyses.

All arthropods were identified to the lowest taxonomic level possible. Five specimens per morphospecies or less were selected for DNA barcoding. One or two legs were removed from each specimen for DNA extraction. Lab work followed standardized protocols for DNA extraction, barcode amplification and sequencing (deWaard et al. 2008). DNA was extracted using a glass-fiber column-based protocol (Ivanova et al. 2006). The primer cocktail C_LepFolF and C_LepFolR (Hernandez Triana et al. 2014) was used to amplify a 658 bp fragment of the cytochrome oxidase I (COI) gene. The PCR thermal regime consisted of an initial denaturation at 94 °C for 1 min; five cycles at 94 °C for 1 min, 45 °C for 1.5 min and 72 °C for 1.5 min; 35 cycles of 94 °C for 1 min, 50 °C for 1.5 min and 72 °C for 1 min followed by a final cycle at 72 °C for 5 min. PCR amplicons were visualized on a 1.2% agarose E-Gel (Invitrogen) and then diluted 1:10 with sterile water. Amplicons (2–5 μL) were bidirectionally sequenced using sequencing primers M13F or M13R (Messing 1993) and the BigDye Terminator v.3.1 Cycle Sequencing Kit (Applied Biosystems, Inc.) on an ABI 3730xl capillary sequencer at the Canadian Centre for DNA Barcoding (CCDB, https://ccdb.ca). All records were uploaded to the Barcode of Life Data System (BOLD, http://www.boldsystems.org) and are publicly available through the dataset DS-BRDNT (doi: dx.doi.org/10.5883/DS-BRDNT) and GenBank (accession: OP587213–OP587227; OP599355–OP599544).

All nests were thawed and then sieved twice using 5 mm and 2 mm sieves respectively (Fig. 1). Nest components larger than 5 mm were manually sorted for arthropod remains; these were collectively stored in 5 ml or 50 ml tubes. Fragments too large to fit through the 2 mm sieve were processed alike. Nest components (i.e. “dust samples”) passing the 2 mm sieve were scooped and placed in 500 ml wide neck Nalgene bottles. Additionally, all feathers found in a nest sample were transferred into a single 5 ml or 50 ml tube. This led to a total of 19 bottles (500 ml), 27 larger tubes (50 ml) and 16 smaller tubes (5 ml) which were stored at -20 °C until further processing in a low-DNA environment (Fig. 2). All work was done on DNA-free surfaces (cleaned with 1.5% bleach and 70% ethanol) and by using DNA-free gloves, which were changed between nests. Sieves and forceps were cleaned with 5% bleach, rinsed with water, and dried after work on a nest.

Arthropod remains obtained during the sieving were ground to fine powder. For larger samples we used 40 ml grinding chambers and an IKA Tube Mill control (IKA, Breisgau, Germany) at 25,000 rpm for 2 × 3 min. Smaller samples were ground in 20 ml homogenization chambers using the IKA ULTRA‐TURRAX Tube Drive Control System (IKA, Staufen im Breisgau, Germany) at 4,000 rpm for 30 min with 10 steel beads (diameter, 5 mm). All samples were covered with a 20:1 mixture of TES (0.1 M TRIS, 10 mM EDTA, 2% sodium dodecyl sulphate; pH 8) and Proteinase K (Qiagen, 20 mg/ml), thoroughly vortexed, and incubated at 56 °C overnight. Subsequently, we used 200 µl of lysate per sample for DNA extraction with the DNeasy kit (Qiagen) following manufacturer’s instructions, except for the elution step, which was done twice using 50 µl AE buffer per sample. Feather samples were lysed with 300 µl of lysis buffer (20: 1 ratio for TES and Proteinase K). Incubation and DNA extraction were identical to the ground arthropod samples.

We generated four (and in exceptions 5) DNA extracts for each nest dust sample: two with the DNeasy PowerSoil Kit (Qiagen) and two with the DNeasy kit followed by processing with the DNeasy PowerClean Pro Cleanup Kit (Qiagen). For four dust samples, three PowerSoil extractions were carried out. First, 2 × 0.25 g dust were processed individually using the DNeasy PowerSoil Kit. Instead of the 10 min vortex step for cell lysis as suggested in the manufacturer instructions, samples were vortexed only briefly and then incubated at 70 °C for 5 min. This process was carried out twice before returning to the remainder of the standard protocol. The remaining nest dust was separated into two 50 ml falcon tubes (in case of large dust samples the amount was limited to 10 g of dust per tube) and covered with a 20:1 mixture of TES and Proteinase K. The lysis and extraction with the DNeasy kit were the same as for the ground arthropod samples, but the generated extracts were additionally subjected to inhibitor removal with the DNeasy PowerClean Pro Cleanup Kit (Fig. 2). For each sample type and extraction, at least one negative extraction control was processed along the samples. Success of the extractions was confirmed by measuring the total DNA concentration with the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific).

All generated DNA extracts were subjected to a two-step PCR protocol (Elbrecht and Steinke 2018) with total PCR volumes of 25 µl using the Multiplex PCR Master Mix Plus (Qiagen), 0.5 mM of each primer, and molecular grade water. The primers BF3 (5’ CCHGAYATRGCHTTYCCHCG 3’) and BR2 (5’ TCDGGRTGNCCRAARAAYCA 3’) and 4 µl of DNA extract were used for the amplification of the COI gene in the first PCR (Elbrecht and Leese 2017; Elbrecht et al. 2019). The second PCR was carried out using fusion primers containing sample-specific inline tags and 2 µl of PCR product generated in the first PCR (Elbrecht and Leese 2017; Elbrecht et al. 2019; Elbrecht and Steinke 2019). The thermocycling conditions for both PCRs were 95 °C for 5 minutes, 35 (1^st PCR) or 24 (2^nd PCR) cycles of 95 °C for 30 seconds, 50 °C for 30 seconds, and 72 °C for 50 seconds with a final extension at 68 °C for 10 minutes. A total of 11 PCR negative controls were used and after PCR 2 the amplification success was checked on 1.5% agarose gels.

For amplicon normalization we processed 20 µl of each PCR 2 product using the SequalPrep Normalization Plate Kit (Invitrogen). We pooled 5 µl per normalized sample into a 5 ml reaction tube. After pooling, its contents were vortexed and distributed into three 1.5 ml reaction tubes. These went through clean up using the SPRIselect Kit (Beckman Coulter) and the Left Side Size Selection procedure with a sample-to-volume ratio of 0.75 and final elution in 30 µl molecular grade water. The eluates of the three reaction tubes were pooled again, the resulting library was checked on an 1.5% agarose gel and its DNA quantity was measured five times with the Qubit dsDNA HS Assay Kit (7.1 ng/µl average). Sequencing was carried out at the University of Guelph’s Advanced Analysis Centre using an Illumina MiSeq with the 600 cycle Reagent Kit v3 (2 × 300) and 5% PhiX spike in. Sequencing results were uploaded to the NCBI Sequence Read Archive (SRA, Genbank, accession: SRR23716567).

Obtained sequences were processed with JAMP (https://github.com/VascoElbrecht/JAMP). During demultiplexing, only sequences with perfectly matching tags were considered for further processing. Forward and reverse reads were merged with U_merge_PE() accessing Usearch (Edgar 2010) requiring 75% of the bases to match. After removing the primer sequences with Cutadapt (default settings; Martin 2011) and a preliminary check of the read length distribution in the remaining dataset, fragments with lengths between 200 and 500 bp were selected for further processing. An expected error value of 2 was used to remove sequences with poor quality as implemented in Usearch (Edgar and Flyvbjerg 2015). During the denoising process, reads were dereplicated per sample and amplicon sequence variants (ASV) with less than 5 reads per sample were removed, as were ASV below 0.001% relative abundance in at least one sample. Additionally, chimeras were identified and removed. The obtained ASV were mapped against a custom sequence database consisting of all publicly available COI sequences on BOLD (as of May 5^th 2021) assigned to “Ontario” using the USEARCH algorithm. All reference sequences were degapped and ambiguous base calls of more than seven Ns in a row were removed prior to mapping. Our primary aim was to identify Annelida, Arthropoda, Chordata, and Mollusca. As the BOLD COI barcode database is very comprehensive for these taxa in southern Ontario, we refrained from separately assigning sequences to different taxonomic levels based on similarity percentages. Only ASV with at least 99% similarity were used and assigned to species level (Elbrecht et al. 2017b). All matched ASV were subjected to a manual plausibility check. This resulted in the use of genus-level information if a) COI was not suitable to distinguish between individual species or b) if the detected species does not occur in the study area, but the genus does. ASV for which no species and/or genus name were available, but which passed this similarity threshold during mapping, were kept in the dataset. For any sequences detected in negative controls, twice the number of reads were removed from all co-extracted samples containing the respective species.

Each of the detected taxa was assigned to a category of reasons for occurring in a bird nest: “bird food” (based on the feeding preferences of the respective nesting species), “feeding in nest” (i.e. saprophagous, herbivore, fungivore), “nest material” (i.e. feathers and fur), “parasite” (i.e. bird ectoparasites and parasites of other taxa) and “visitor” (e.g. mice or squirrel visiting the nest). Additionally, each taxon was associated with a consumer type based on the behavior at its adult life stage: “decomposer”, “not eating”, “palynivore”, “fungivore”, “herbivore”, “omnivore”, “parasite”, “predator”, “saprophagous” and “scavenger” (Hicks 1959, 1962, 1971; Levesque-Beaudin et al. 2020; https://Bugguide.net). For both categorisations, an unambiguous assignment was not always possible.

All calculations and visualizations were made in R Version 4.2.2 (R Core Team 2021) using the packages “ggplot2” (Wickham 2016), “dplyr” (Wickham et al. 2022), “ggrepel” (Slowikowski 2021) and “ggpubr” (Kassambara 2020). Poisson regressions with a log-link function were used to test for significant differences in the number of detected taxa between the four different extract types (ground arthropods, dust DNeasy extracted, dust PowerSoil extracted, feathers). The unbalanced distribution of sampled nests between bird species did not enable the application of ordination methods to examine differences between species communities with respect to the nesting bird and resulted in the use of descriptive statistics for the obtained data.

The Miseq run produced 27,407,037 reads and 9.69% of these were discarded during demultiplexing. Of the discarded sequences, 19.4% were assigned to PhiX using a similarity threshold of 90%. On average, 76.5% of the forward and reverse reads per sample could be merged and the removal of the primer sequences was successful for 99% of the reads. After length selection, primarily removing reads of primer-dimers, 2,523,299 reads, on average 14.5% per sample, remained for further processing. During the denoising 1,393,893 reads were removed (including chimeras) and 791,123 reads were used for taxonomic assignment. The extraction controls processed alongside the ground arthropods and the feather samples showed contaminations with Turdus migratorius, Lumbricus terrestris Linnaeus, 1758, Lymantria dispar Linnaeus, 1758, and Ornithonyssus sylviarum Canestrini & Fanzago, 1877. After down-correcting the reads in the affected samples, removing Homo sapiens Linnaeus, 1758 sequences, and a general plausibility check, 352,692 reads belonging to 103 taxa were used for the analysis.

Ground arthropod samples and PowerSoil nest dust extractions each accounted for 39% of the reads, whilst DNeasy nest dust extractions and feather samples accounted for 13% and 9%, respectively. Aves (68%) and Insecta (23%) were the two classes to which most reads were assigned to (Fig. 3). Of the 103 taxa detected with the metabarcoding approach, 11 taxonomic assignments were to family level, 10 to genus level, and 82 to species level, respectively. Hereafter, the word “taxon/taxa” is used to describe distinct taxonomic entities in an effort to include the cases in which reads were assigned to a family or genus but could not be matched to any of the species included in the database.

Figure 3.

The bubble plot shows the proportion of reads per bird nest which are associated to the four different sample/extraction types (“a”: ground arthropods, “dd”: DNeasy extraction, “dp”: PowerSoil extraction, “f”: feather samples), i.e., proportions for each nest sum up to 1. Nests are sorted and bubbles are coloured according to the nesting bird species and the proportional reads are separately displayed for all detected taxonomic classes.

For nine nests, taxonomic assignments were possible for small and medium ground arthropod subsamples picked from the 5 mm and 2 mm sieve, respectively. Of the taxa detected in these samples (mean 9.11 ± 5.86 SD), an average of 5.5% were simultaneously detected in both subsamples. Both subsamples of the DNeasy extracted nest dust led to taxonomic assignments for 10 nests, with a mean of 3.10 ± 2.02 SD taxa of which on average 75.50% were detected in both subsamples. The PowerSoil extraction of nest dust enabled taxonomic identifications in more than one subsample for 17 nests. Of the detected taxa (mean: 4.82 ± 2.24 SD), 47.31% were on average detected in more than one subsample per nest.

The analysis of nest dust resulted in 19 nests with taxonomic assignments in at least one subsample based on the PowerSoil extraction, 13 had DNeasy-based detections, and 62.86% of these (DNeasy detections) were observed with both methods. Per nest, significantly more taxa were detected from dust with the PowerSoil extraction method compared to the DNeasy extraction method (Table 2; mean PowerSoil: 4.89 ± 2.28 SD; mean DNeasy: 2.77 ± 2.05 SD). A comparison of the number of detected taxa per nest and sample type revealed significantly less detected taxa in DNeasy extractions compared to the other three sample types (i.e. PowerSoil extracted nest dust, ground arthropods and feathers; Table 2; Fig. 4).

Figure 4.

Box plots of the detected taxa (lowest possible taxonomic level) in the four sample types: ground arthropods (a), DNeasy extracted nest dust (dd), PowerSoil extracted nest dust (dp), and feather samples (f). Significantly less taxa were detected after DNeasy extraction of dust (Table 2).

In the four nests processed using emergence traps, 24 distinct taxa (9 different orders) were detected of which 58% were identifiable to species level using DNA barcoding. Insecta made up 92% of the detected taxonomic diversity in emergence traps with Diptera, Coleoptera and Psocodea each contributing 25%, 25%, and 21%, respectively (Fig. 5). The metabarcoding analysis of ground arthropods, nest dust, and feather samples led to the detection of 103 distinct taxa, 82 of which were identifiable to the species level (See Supporting Information 1 for the full data table). Arthropoda accounted for 72% of the distinct taxa and at the order-level, Diptera (class Insecta) and Passeriformes (class Aves) accounted for the highest percentage of taxonomic diversity (18% and 14%, respectively, Fig. 5). Regarding the arthropod diversity detected via metabarcoding, 53 taxa (38 species, 5 genus-level identifications, 10 family-level identification) only occurred in one nest. These unique detections constituted 23.5% of the total arthropod reads and were primarily (79%) Insecta. Most arthropod reads (20.9%) were assigned to Protocalliphora sp. Hough, 1899 which was detected in 7 nests. These bird blowflies lay their eggs in bird nests and parasitize nestling (Sabrosky et al. 2011). Of the 21 taxa (15 species, 5 genus-level identifications, 1 family-level identification) detected in more than one nest, Niditinea orleansella Chambers, 1873 and Ceratophyllus sp. Curtis,1832 occurred in 10 nests each. The moth genus Niditinea Petersen, 1957 is known to occur in bird nests feeding on nest material and dried animal matter (Townsend and Waring 2018), and Ceratophyllus is a genus of fleas, for which several species are known to be associated with birds (Tripet and Richner 1997; Brinck-Lindroth and Smit 2007).

Figure 5.

Taxonomic diversity (i.e., proportion of distinct taxa) detected from a) the metabarcoding of nest materials from 20 nests b) the nests placed in emergence traps (n = 4) and c) the metabarcoding of nest materials previously placed in emergence traps (n = 4; only arthropod detections displayed). The metabarcoding data in panel a are based on 121 extracts.

For the four nests which were placed in emergence traps and subsequently analyzed via metabarcoding, the mean arthropod richness detected via metabarcoding was 12.50 ± 7.94 SD taxa compared to 8.00 ± 3.16 SD taxa calculated from emergence trap detections. The most frequently detected order was Diptera (class Insecta) in the metabarcoding data and Coleoptera (class Insecta) in the emergence trap data (Fig. 5). Of the 40 taxa detected via metabarcoding, 52.5% were classified as “bird food” or ambiguously “feeding in the nest or bird food”, indicating that these taxa would unlikely be detected in an emergence trap. Five taxa (Niditinea, Ceratophyllus, Protocalliphora, Gaurax pallidipes Malloch,1915, and Liposcelis corrodens Heymons, 1909) were detected both in emergence traps and with DNA metabarcoding. Niditinea, Ceratophyllus, and Protocalliphora were the most frequently detected taxa in all 20 nests, and both G. pallidipes and L. corrodens were detected more than once via metabarcoding. Altogether, 25 bird species, all of which occur in the study region, were detected with metabarcoding. Except for the Black-Capped Chickadee (Poecile atricapillus Linnaeus, 1766), the nesting bird species was always detected in the respective nest(s) and accounted for 32–71% of the total Aves reads (summed up per nesting species; Fig. 6). Most bird species (n = 18) were detected in the eight House Wren (Troglodytes aedon Vieillot, 1809) nests.

Figure 6.

The bubble plot shows the proportion of reads which are assigned to one of the 25 detected bird species for every nesting bird species i.e. proportions for each nesting bird sum up to 1. Detections of the nesting bird species itself are coloured in orange, detections of species which are foreign to the nest, are coloured in blue.

Of the 103 taxa detected with DNA metabarcoding, 73% could be associated with a specific reason for occurring in the nest samples: the occurring Arthropoda were primarily food for the nesting birds or feeding in the bird nest (Fig. 7). With respect to consumer type, 68% of the detected taxa had an unambiguous link at the adult stage and most of the detected arthropod taxa could be associated with herbivory (Fig. 7).

Figure 7.

Bar plot showing the number of arthropod taxa which could be unambiguously assigned to a) a reason for occurring in a bird nest (n = 74) and B) a consumer type (n = 69).

The authors have declared that no competing interests exist.

No ethical statement was reported.

Laboratory work was supported by the Canada First Research Excellence Fund and represents a contribution to the University of Guelph’s “Food from Thought” program.

VLB conceived the original idea, which was developed further together with DS and BT. DS secured the necessary funding. VLB collected the bird nests and was responsible for the emergence traps and the consecutive morphological identification of emerging arthropods. VLB and BT sieved the nests and sorted the arthropod remains, BT was responsible for the consecutive laboratory processing of the generated samples which was carried out together with MB. BT was responsible for analyzing the generated data, VLB wrote a first draft of the manuscript to which all co-authors contributed critically before final approval for publication.

Valerie Levesque-Beaudin https://orcid.org/0000-0002-6053-0949

Dirk Steinke https://orcid.org/0000-0002-8992-575X

Bettina Thalinger https://orcid.org/0000-0001-9315-8648

The raw metabarcoding data has been uploaded to GenBank SRA (accession: SRR23716567). All data obtained from emergence traps and molecular data used for the consecutive analyses have been uploaded to Figshare and are available at https://doi.org/10.6084/m9.figshare.22203616.v1.