Bitter fruits of hard labour: diet metabarcoding and telemetry reveal that urban songbirds travel further for lower-quality food

Rapidly increasing urbanisation requires mitigation against associated losses of biodiversity and species abundance. In urban-breeding birds, altered food availability for nestlings is thought to reduce reproductive success compared to forest populations. To compensate for shortages of preferred foods, urban parents could increase their search effort for optimal diets or provision other foods. Here, we used telemetry and faecal metabarcoding on blue tits from one urban and one forest populations to compare parental effort and comprehensively describe nestling diet. Urban parents travelled on average 30% further than those in the forest, likely to offset limited availability of high-quality nestling food (i.e. caterpillars) in cities. Metabarcoding, based on a mean number of 30 identified taxa per faeces, revealed that the diets of urban chicks were nonetheless substantially shifted to include alternative foods. While in the forest caterpillars comprised 82 ± 11% of taxa provisioned to nestlings, in the city they constituted just 44 ± 10%. Pre-fledging chick mass as well as offspring numbers were lower in urban than in forest-reared broods. Thus, at least in our comparison of two sites, the hard labour of urban parents did not fully pay off, suggesting that improved habitat management is required to support urban-breeding birds. Electronic supplementary material The online version of this article (10.1007/s00442-020-04678-w) contains supplementary material, which is available to authorized users.

of hatchlings/eggs, fledging success as the maximal numbers of fledglings/hatchlings, and used fledging mass (mass on post-hatching day 13) as a proxy for the probability of survival. After fledging, we recorded any dead chicks.
At each site, 8 focal nestboxes were chosen according to their suitability for telemetry, the ideal terrain being as even as possible, with a slight slope and some potential vantage points. One brood in the city died at day 7 of nestlings' lives; for this brood we did not collect nestling mass data, faecal samples or video footage. We caught one of the parents from each nestbox on the 5th (n=9), 6th (n=6) or 7th (n=1) day of the nestlings' lives while it provisioned its brood. The parents were tagged with a radiotransmitter (see below), and weighed, ringed, sexed and aged.
Infrared camera systems were installed at each nestbox on day 7 and 11 of the nestlings' life (forest n=8, city n=8; custom-made systems by the University of Glasgow's Bioelectronics Unit;Pollock et al. 2017). On the 13th day of the nestlings' lives, we collected at least two faecal samples per nest. Samples were collected by picking up the chick and immediately holding a vial to its cloaca, in some occasions gently massaging the lower abdomen. Faecal samples were collected in vials containing 1.25ml of 100% ethanol and stored at -20°C during the field-season. Chicks were then weighed, ringed, and a blood-sample taken.
Telemetry: We tagged adults with radiotransmitters as described in Nord et al. 2016, using PIP31 single celled tag (Biotrack, Dorset, UK; weight 0.35g; maximal dimensions 7x7x4mm; battery type Ag337). Briefly, we trimmed feathers on the birds' lower back and attached a radiotransmitter using eyelash glue (Duo Eyelash Adhesive Clear White, waterproof) and a small amount of superglue. The birds were tracked over the next 1-4 days with Lotek wireless SRX400 receivers and Yagi antennas by two observers (CJ and HM). To triangulate the position of the bird, observers stood at an approximate 90⁰ angle from the estimated bird location and took compass bearings in the estimated direction of the bird every 2 min, for periods of 30 min. The locations where the observers stood were recorded using GPS (Garmin, eTrex 30x). Whenever the bird was seen at any of the 2-min points during the tracking period, a bearing and an estimated distance from the observer were recorded. Observers coordinated their bearings with walkie-talkies, and synched watches. We carried out a total of 3-5 tracking periods of 30 min per bird, recording a total of 666 position fixings spread out across the day.
For quality control of the data, the quality of the signal was recorded for each observation, based on the volume and stability of the signal. Categories used to indicate quality were "good" (when there was a steady volume signal with a clear direction), "moving" (when direction was detectable but volume varied unexpectedly), and "bad" (when signal was faint and direction hard to detect). Any fixes where the signal quality was noted as being "bad" were excluded from analysis; there were more "bad" fixings in the city than the forest (45 and 26 respectively), likely due to interference with buildings. Triangulation calculations were carried out using the Sigloc package (Berg 2015) within R 3.3.1. Once all bird positions had been calculated, points were plotted using Google Earth (version 7.0.2.8415, 2012) and points deemed unlikely based on visual assessment of the bird's position on the map were eliminated. The criteria used for exclusion were: if the point was a clear outlier compared to the majority of the data, or if the actual location of the bird was considered unrealistic (such as over water or behind rows of buildings). The number of position fixes after cleaning up the data was 570 (city: n=303; forest: n=267). We used the package Geosphere (Hijmans et al. 2012) to calculate foraging distances, considered as the distance (m) between the nestbox and each location of the bird, in a straight line.
Video recording of parental provisioning: Infrared camera systems were focused on the entrance hole from inside the nestbox and recorded for 24 hours, starting at approx. 5pm. Four standardised half-hour periods (8:00-8:30 and 19:00-19:30 per sampling day) were extracted for each nestbox using VideoLAN VLC Media Player. From counted parental entries, we calculated provisioning rate as the number of parental entries into the nestbox per half-hour period.
We classified provisioned prey items into two categories: caterpillars and non-lepidopteran invertebrates. Where possible we identified invertebrates to family level (e.g. aphids). We were unable to identify 16% (0.16±0.17) of all items delivered, with no bias for site (p>0.05). These items were discarded for analysis. For each half-hour period of footage, the number of items from each category delivered was divided by the total number of visits to obtain proportions of visits. The dimensions of the caterpillars delivered was estimated by obtaining total length (L) and mean width (W) using the diameter of the nest hole as a reference, and the formula (π/4)*L*W 2 (Blondel et al. 1991) to calculate the volume.

Metabarcoding and bioinformatics
PCR. PCRs were conducted in 15μL reactions containing 1X Multiplex PCR Mastermix (Qiagen), 0.13μM of each primer, and 3uL DNA extract. PCR thermocycling conditions followed an initial denaturation of 95°C for 5 min; 35 cycles of 95°C for 30 sec, 50°C for 1 min, and 72°C for 30 sec; and a final extension of 68°C for 10 min. For rbcL, multiplexed PCR amplification was conducted by simultaneously using all three primer pairs with the same mastermix and thermal cycling protocols as described above. PCR was conducted up to seven times for each locus, with results examined on 2% agarose gels stained with ethidium bromide, until three positive PCR products were obtained.
Library construction and validation. After cleaning and pooling the PCR products as indicated in the main text, a dual indexing PCR was performed: 50μL reactions contained 1x HIFI Master Mix (KAPA), 0.5μM of each indexed primer, and 10uL bead cleaned PCR product. The thermal cycling conditions were: 72°C for 3 min; 98°C for 30 sec; 12 cycles of 98°C for 10 sec, 63°C for 30 sec, followed by 72°C for 3 min. Reactions were cleaned, quantified using the Qubit High Sensitivity DNA assay, and then pooled in equimolar ratios. The sequencing pool was analysed on a 2100 Bioanalyzer (Agilent Technologies) using the high sensitivity DNA kit. The final pool concentration was determined via quantitative PCR with the KAPA Library Quantitative Kit.
Bioinformatics. Unless otherwise specified, default parameters were employed. Raw sequences were analysed using FastQC (Andrews et al. 2011) and adapters were trimmed using CutAdapt v1.10 (Martin 2011). Following Schirmer et al. (2015, quality trimming was conducted using Sickle v1.33 (Joshi & Fass 2011) with the -x option to prevent 5' trimming (Joshi & Fass 2011). Then error correction was performed using BayesHammer (Nikolenko et al. 2013) through the SPAdes program v3.10.1 (Bankevich et al. 2012), and reads were merged using PEAR v0.9.6. (Zhang et al. 2014). Data for each primer set were split using a custom python script and PCR primers were trimmed off using CutAdapt. For the COI dataset, we filtered non-target sequences (e.g. those potentially belonging to the birds or humans) using BLAST: first, identical reads for each sample were collapsed using the obiuniq function in obitools (Boyer et al. 2016), while retaining the read counts. They were then queried against the curated COI Midori-LONGEST database (Machida et al. 2017) using blastn in BLAST+ v2.6.0. A custom python script was used to retain reads, and associated abundance information, for which the best match (based on bitscore) was classified as belonging to phylum Arthropoda, the alignment was ≥50% the length of the query sequence, and the e-value was ≤ 10 -5 . The data were then filtered for chimeric sequences using usearch v11.0.667 (Edgar 2010) implemented through a modified version of the program DAMe (Zepeda-Menoza et al. 2016) available at https://github.com/shyamsg/DAMe. The DAMe script convertToUSearch.py was used prepare for OTU clustering, and to filter sequences ± 30bp of the target sequence length (without primers).
OTU clustering and taxonomy assignment. As discussed by Clare et al. (2011) and outlined by Razgour et al. (2011), we tested OTU clustering similarity thresholds from 94 -100%, and used plots to identify the best threshold (97%) where over-and under-splitting were minimized (Krüger et al. 2014;Trevelline et al. 2016). We used the tabulateSumaclust.py script from the DAMe package to make OTU tables, retaining only OTUs with ≥ 5 sequences. To assign taxonomy, we used a BLAST search of the Genbank NT database: Using blastn in BLAST+ 2.7.1 we retrieved the best 20 matches for each representative OTU sequence. Using custom python and R scripts (Alberdi et al. 2018;Aizpurua et al. 2018), we retained the highest common taxonomic information among matches with the maximum bitscore, and assigned taxonomy to each OTU based on identity: For matches with ≥95% identity we assigned orderlevel taxonomy; ≥96.5% we assigned family-level, and for ≥98% we assigned genus and species-level taxonomy.