Metabarcoding-based dietary analysis of hen harrier (Circus cyaneus) in Great Britain using buccal swabs from chicks

The hen harrier is a heavily persecuted bird of prey in Great Britain since its diet includes Red grouse, a game bird shot in driven and walked-up grouse shooting. Unlike walked-up shooting where shooters walk up moors and flush grouse with dogs, in driven grouse shooting the grouse are driven by beaters towards static shooters. Driven grouse moors (DGMs) are increasingly being managed to sustain high densities of Red grouse intensifying a long-standing conservation conflict between conservationists and grouse moor keepers. A metabarcoding approach was used with degenerate universal cytochrome oxidase I and cytochrome b primers along with hen harrier blocking primers. A novel sampling method was used to detect prey in buccal swab samples from chicks from broods across Great Britain from both managed and unmanaged moorland habitats. This resulted in detection of 62 species of prey across 51 broods with Meadow pipit, Red grouse, Wren, Skylark, and voles being most frequently detected. Frequency of occurrence data and species accumulation curves reveal high incidence of Red grouse and low prey species richness in the diet of hen harriers in DGMs but low incidence of Red grouse and high prey species richness in walked-up and unmanaged moors. Waders were only detected within walked-up and unmanaged moors and not within DGMs where they have been reported to occur at high densities. Regional species detected included endemic species such as the Orkney vole seen only in Orkney. This study represents the first metabarcoding-based dietary analysis in a raptor using buccal swabs.


Introduction
The hen harrier (Circus cyaneus) is a raptor species with a wide distribution in the northern hemisphere. The IUCN has assessed its conservation status as 'least concern' although there is evidence of a global population decline due to habitat loss as a result of intensified agriculture, disappearances of marshes, and afforestation programs (BirdLife International 2016). Fielding et al. (2011) estimated that there was enough habitat to sustain a population of 2036-2390 breeding pairs of hen harriers in Great Britain (1467-1790 in Scotland, 323-340 in England, and246-260 in Wales), but this has never been reached since severe persecution especially in northern England and southern and eastern Scotland, has consistently limited numbers (Etheridge et al. 1997;Sim et al. 2007). The latest population survey conducted in 2016 (Wotton et al. 2018) showed an estimated 545 pairs in the whole of the United Kingdom (excluding the Isle Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1059 2-019-01215 -y) contains supplementary material, which is available to authorized users.
Hen harriers are persecuted in Great Britain since their diet includes the game bird, the Red grouse (Lagopus lagopus scotica). Red grouse are shot in large numbers in driven grouse shooting during the shooting season, with an average of 72 brace (144 birds) per day of shooting (Sotherton et al. 2009). This type of grouse shooting where grouse are driven by beaters towards the shooters requires the grouse population to be far higher than would naturally occur e.g. 5.6-7.74 willow grouse (Lagopus lagopus) per km 2 have been reported in unmanaged habitat in Sweden by Willebrand et al. (2011) and 104.9 Red grouse per km 2 have been reported in driven grouse moors (DGMs) in northern England by Robertson et al. (2017). Dense populations are maintained through land management which involves heather burning and predator control (Thirgood et al. 2000). Predator control limits numbers of naturally occurring predators not protected by legislation e.g. stoats (Mustela erminea), weasels (Mustela nivalis), foxes (Vulpes vulpes) and crows (Corvus spp.). Management has been intensified on some DGMs in the last few decades with burning rotations getting shorter, grouse being medicated against a gastrointestinal nematode using anthelmintics in grit, and non-predatory species which could indirectly harm grouse being culled e.g. mountain hares which carry ticks that transmit Louping ill virus to Red grouse (Thompson et al. 2016). An alternative form of grouse shooting, walked-up shooting, involves shooters walking up moors and shooting grouse that are flushed using dogs. This form of grouse shooting, practised in Wales, is less intensive in its management practices and results in far smaller numbers of grouse being shot (Sotherton et al. 2009).
Previous dietary studies have shown that Red grouse can be a major part of the diet of hen harriers on DGMs (Redpath and Thirgood 1997). This has resulted in land managers and gamekeepers perceiving any degree of predation as an economic loss to the shooting estate (Hanley et al. 2010). Despite being listed on the Red List of Birds of Conservation Concern and protected by The Wildlife and Countryside Act of 1981, hen harriers continue to be affected by persecution. Using satellite tracking data from 58 hen harriers, Murgatroyd et al. (2019) reported a high likelihood of mortality with increased use of grouse moors.
Grouse moor management has been reported to improve breeding success and abundance of ground nesting species other than Red grouse e.g. Tharme et al. (2001) reported higher breeding densities of waders such as golden plover (Pluvialis apricaria), lapwing (Vanellus vanellus), and curlew (Numenius arquata) on managed moors compared to unmanaged moors. Further, Fletcher et al. (2010) reported a threefold increase in breeding success of golden plover, lapwing, curlew, Red grouse and meadow pipit (Anthus pratensis) on moorland as a result of predator control measures reducing numbers of foxes and carrion crows (Corvus corone), and Ludwig et al. (2019) reported similar results where a positive correlation was observed between predator control as part of grouse moor management and numbers of some waders. However, a study on the relationships between wader declines and upland land use by Amar et al. (2011) reported that population declines of Golden Plovers were greater in moors intensively managed for grouse shooting. These authors found no evidence for declines in Curlew populations being lower in managed moorlands. Concern has also been growing in recent years over the negative environmental impacts of intensive management practices used on grouse moors. Frequent illegal killing of protected raptors on DGMs, negative environmental impacts of frequent burning and anthelmintic usage, and a likely increase in incidence of respiratory cryptosporidiosis in Red grouse due to high population densities in DGMs are all cited by Thompson et al. (2016) as reasons to question the sustainability of driven grouse shooting. An alternative view was however, presented by Sotherton et al. (2017) who whilst acknowledging that illegal killing of raptors cannot continue, highlighted the positive effects of managed burning on protecting dwarf shrub heath habitat and prevention of wildfires. They also questioned the perceived environmental impact of anthelmintic drugs in grit. Thus, opinions remain divided and the conservation conflict continues (Hodgson et al. 2018).
There are many examples of dietary analysis playing an important role in the conservation of species (Chase Grey et al. 2017;Aziz et al. 2017;Buglione et al. 2018). In case of the hen harrier, although the diet has been relatively well studied in mainland Scotland and Orkney, most of these earlier studies predate the intensification of grouse moor management in mainland Scotland (Balfour and Macdonald 1970;Picozzi 1978Picozzi , 1980Redpath 1992;Redpath and Thirgood 1997;Redpath et al. 2001). Observational studies and morphological analysis of pellets showed that the diet of hen harriers includes passerines, game birds, Lagomorpha, and small mammals such as the short-tailed field vole (Microtus agrestis). Meadow pipit and Red grouse were found to be the most abundant species in mainland Scotland with 23-45% and up to 34% of the prey items brought to the nest respectively (Picozzi 1978;Redpath and Thirgood 1997). Only up to 7% of the prey brought to the nest were identified as Red grouse in Orkney where grouse densities are low (Balfour and Macdonald 1970;Picozzi 1980). Redpath (1992) found a positive correlation between red grouse density and prevalence in hen harrier diet i.e. if grouse populations were low in density, hunting shifted to passerines. These earlier diet studies may not necessarily represent current availability of prey species, making it important to reinvestigate the hen harrier's diet across different regions.
Sensitive DNA based techniques provide the means to investigate diet across a wide geographical range whereas conventional methods are labour intensive, range restricted, and could contain observer biases when species are wrongly identified or absent. Using DNA has repeatedly been shown to be a reliable method in identifying dietary species from stomach and faecal samples using either prey specific primer assays (e.g. Thalinger et al. 2016), universal primers in combination with cloning and Sanger sequencing (e.g. , or more recently, High-Throughput Sequencing (HTS) (e.g. Jarman et al. 2013). HTS has greatly improved sensitivity and streamlined the analysis of complex prey mixtures using universal primers. Use of blocking oligonucleotides containing C3-spacers overlapping the 3′ end of one of the two primers, preventing amplification of the target region in the predator species (Vestheim and Jarman 2008), also enables blocking of predator DNA amplification when using universal primers, which is especially useful in diet analysis of predators that feed on closely related species. We used a novel approach for dietary analysis in birds by using DNA from buccal swabs from hen harrier chicks. Samples were amplified using two universal primer sets designed to amplify short fragments of the COI and Cytb genes across all terrestrial mammals, birds, and reptiles resident/breeding in the UK. Using HTS based metabarcoding, we investigated the diet of hen harriers from across several different regions in Great Britain within both managed and unmanaged moorland habitats. Our specific aims were to: (1) assess how suitable buccal swabs are for detection of prey species; (2) use the Frequency of Occurrence (FOO) i.e. the proportion of samples in which a particular prey species is detected, as a measure of differences in prey species in the diet of hen harriers across different moorland habitats; (3) confirm higher FOO for Red grouse in DGMs compared to unmanaged moorland habitats; and (4) use species accumulation curves as overall indicators of prey species richness in the diet of hen harriers.

Materials and methods
Samples DNA samples extracted from buccal swabs were available from a previous study (van Hoppe et al. 2016) and a subset was chosen for dietary analysis. All samples were collected from 22 to 30 day old hen harrier chicks by handlers with the relevant Schedule 1 license during 2006, 2007, 2008, and 2014 from England (20 broods, 5 sites, mean broods per site 3.3), mainland Scotland (14 broods, 9 sites, mean broods per site 1.5), Orkney (12 broods, 2 sites, mean broods per site 6), and Wales (7 broods, 3 sites, mean broods per site 2.3). Samples were obtained from broods within/< 3 km from DGMs, within walked-up moors and from unmanaged moors. Two DGMs from mainland Scotland were no longer managed for driven grouse shooting at the time of sampling and were therefore, classified as historical DGMs. Samples originated from 51 different broods with up to five individuals per nest. One individual was selected from each brood (sample with the highest amount of DNA) to get a broad overview across locations but from four broods (two in England, one in mainland Scotland, one in Orkney), every individual was used.

Controls
Sequencing was performed using three Ion Torrent semiconductor sequencing chips: two 314™ Chips v2 (Ion Tor-rent™) (Chips 1 and 2) and one 318™ Chip v2 (Ion Tor-rent™) (Chip 3). Positive control samples were created to run alongside samples. They contained DNA from hen harrier and five other species available in the laboratory and chosen based on the melting temperature (Tm) of degenerate primers (see Table S1 for details). Species unlikely to be prey such as grey squirrel (Sciurus carolinensis) and hedgehog (Erinaceus europaeus) were added because they were available in-house and displayed suitable Tms to test. The ratio of hen harrier to prey was 10:1 for chip 1 (final concentration 3 ng/µL hen harrier and 0.3 ng/µL total prey) and 20:1 for chips 2 and 3 (final concentration 3 ng/µL hen harrier and 0.15 ng/µL total prey). DNA was quantified using gel electrophoresis alongside Lambda DNA dilutions (chip 1) or Qubit™ dsDNA BR Assay Kits (ThermoFisher Scientific, Waltham, USA) (chips 2 and 3) for addition in equimolar concentrations to the library. Positive control samples were created on separate days to amplicon library construction to prevent contamination. Two extraction negatives generated during hen harrier buccal swab extractions in 2016 were used (one for chip 1 and one for chips 2 and 3).

Library construction and sequencing
Universal degenerate primers were designed for COI and Cytb (see section on primer design in supporting information and Table S2). Sequencing libraries were constructed using the fusion primers approach for the Ion Torrent sequencing platform. The forward fusion primers contained the A adaptor sequence followed by a validated IonXpress barcode, barcode key (GAT), and the universal degenerate forward primer designed in this study. The reverse primer contained the Trp1 adaptor sequence, IonXpress barcode, followed by the universal degenerate reverse primer designed in this study. Five unique barcodes were used in both the forward and reverse primer which allowed sequencing of 25 samples per target on a single chip (Clare et al. 2014) (Table S3). Each sample was amplified with both COI and Cytb primers in separate reactions.
PCRs contained 1X platinum multiplex PCR master mix (Applied Biosystems™, Foster City, CA, USA), 1 µM of the forward fusion primers, 1 µM of the reverse fusion primers, 1 µM of blocking primer, and approximately 3 ng of buccal swab DNA in 12.5 µL reaction volumes. PCR was set up in a flow cabinet and all equipment was UV cross linked for 20 min prior to PCR setup. The following thermocycling conditions were used; 95 °C for 2 min, 5 cycles of 95 °C 30 s, 51 °C for 90 s, and 72 °C for 60 s, followed by 30 cycles of 95 °C for 30 s, 47 °C for 90 s, and 72 °C for 60 s, with a final extension at 72 °C for 10 min. Each chip contained 22 samples, one extraction negative, one PCR negative, one positive control without blocking primers, and one positive control with blocking primers. All PCR products were first visualised by loading 2.5 µL of the PCR product combined with 2.5 µL water and 1 µL 6X tracker dye on a 1.5% agarose gel.
Libraries for chips 1 and 2 were quantified using serial dilutions of Lambda DNA on an agarose gel, diluted to the lowest concentration of PCR product available, and then pooled in equimolar concentrations by using 1 µL of each diluted PCR product. When no amplicons were seen on the agarose gel and in the case of extraction and PCR negatives, 2 µL (volume of the least diluted sample) was added to the pool. Pooled samples were purified using MicroClean (Microzone Ltd., Haywards Heath, UK) according to manufacturer's instructions. Libraries for chip 3 were individually purified using the GenCatch™ Advanced PCR Extraction kit (Epoch Life Science, Sugar Land, TX, USA) according to the manufacturer's instructions, visualised on an agarose gel and quantified using either the Qubit™ dsDNA HS Assay Kit (ThermoFisher Scientific, Waltham, USA) or Qubit™ dsDNA BR Assay Kits (ThermoFisher Scientific, Waltham, USA) depending on band brightness. DNA was then diluted and pooled. Libraries for chip 1 were diluted to 100 pM prior to emulsion PCR but to 200 pM for chips 2 and 3 to compensate for primer dimers and non-specific amplicons still present after purification.
Ion PGM™ Hi-Q™ OT2 kit (Ion Torrent™) was used for emulsion PCR on the Ion OneTouch according to manufacturer's recommendations and 10 µL Ion PGM™ Calibration Standard (Ion Torrent™) added during the preparation of the amplification solution. Sequencing was performed with the Ion PGM™ Hi Q™ Sequencing Kit (Ion Torrent™) according to manufacturer's instructions on an ION Torrent PGM (ThermoFisher Scientific, Waltham, USA).

Sequence reference database
Local reference sequence databases were created for COI and Cytb for species identification by downloading all available potential prey sequences from BOLD and Gen-Bank respectively. Only the target length was used and all identical sequences from the same species were removed. Makeblastdb command was used in the Linux BLAST+ 2.6.0 package (Camacho et al. 2009) to convert the created FASTA files to a BLAST database. The COI database comprised 138 species (including human) and 438 sequences. The Cytb reference database consisted of 604 sequences representing 136 species. Sequences unavailable in BOLD and GenBank were added from samples in house where available (Table S1). The target length for both COI (130 bp including primers) and Cytb (123 bp including primers) showed enough taxonomic resolution to identify all Operational Taxonomic Units (OTUs) to a single species in the database based on more than 97% homology threshold, except in the following cases: two species of Loxia (Loxia curviostra and L. pytyopsittacus) had identical COI and Cytb sequences; two species of crow (Corvus corne and C. cornix) displayed > 98% homology with both markers and three species of Wagtails (Motacilla alba, M. flava, and M. cinerea) displayed homology ranging between 99 and 98% between species.

Species assignment and sequence filtering
BAM files containing IonXpress barcodes were converted to FASTA files using the 'Convert, Merge, Randomise' tool on the online Galaxy platform (Goecks et al. 2010). FASTA files were re-formatted to single line sequences using the 'FASTA Width' tool. The 'Barcode Splitter' tool (allowing for 2 mismatches) was then used to split the reads based on the barcode in the reverse fusion primer. These files were further analysed with the Linux package OBITools v1.2.11 (Boyer et al. 2016). ngsfilter was used to remove primer sequences and the barcode sequence in the reverse primer. Obiuniq was used to cluster identical sequences. Redundant information in sequence headers was removed with obiannotate so that only the read identifier and the read count remained. Obiclean was used to compensate for PCR and sequencing errors, whereby reads of the Obiuniq clusters which were > 97% similar were added together. Obigrep was used to discard sequences with frequencies below 2 reads and sequences < 85 or > 95 nucleotides for COI sequences, and < 72 or > 82 for Cytb to eliminate non-target length sequences (target length 91 bp for COI and 76 for Cytb).
Clusters were assigned to species by comparing filtered fasta files against the custom-made databases. Reads were identified to species when there was > 97% homology to a reference sequence and more than 3% divergence to the next highest matching species with one exception: Meadow pipit and Rock pipit (Anthus petrosus) differed by only 1.1% with COI but > 3% with Cytb. To make sure Meadow pipit and Rock pipit reads were not clustered together during the obiclean step, reads were checked prior to obiclean to make sure prey clusters contained only one species. Reads with identity levels < 97% and the clusters which could not be matched to any species in the local database were run through NCBI's BLAST. When a highest match of < 97% was observed, resulting BLAST alignments were inspected to see whether this was caused by deletions within mononucleotide stretches (common with low quality sequences), and if so, reads were added to the OTU already assigned in the sample.
Prey reads detected in the PCR negative controls were used to subtract reads of those species on that chip. In addition, all clusters identified to species with a relative abundance below 0.5% (calculated using all reads assigned to species) were removed to account for differences in read coverage between chips. All clusters at > 0.5% relative abundance which matched to the following species were removed since they were commonly handled in the laboratory: Pig (Sus scrofa), Sheep (Ovis aries), dog (Canis lupus), Rainbow trout (Oncorhynchus mykiss), Sterlet sturgeon (Acipenser ruthenus) and European sturgeon (Huso huso), red deer (Cervus elaphus), roe deer (Capreolus capreolus), Cow (Bos taurus), and Human (Homo sapiens).

Data analysis
Prey reads were converted to binary presence/absence data. Species accumulation curves were created for all samples with COI, Cytb, and the combined dataset using the specaccum function Vegan (Oksanen et al. 2019); Exact method, 999 permutations in R 3.5.1 (R Core Team 2018). In broods with multiple samples, the sample with the highest amount of DNA was chosen. The broods were categorised according to region (England, Mainland Scotland, Orkney, and Wales) and moor type (DGM, Unmanaged, and Walked up). For moor-type comparisons and Frequency of Occurrence (FOO) calculations, broods containing less than 400 reads identified to species were removed since at least 400 reads are necessary to obtain OTUs with 2 reads (i.e. 0.5% relative abundance) and one individual from each brood was selected as above. FOO was calculated for samples grouped by moor type (number of samples detected/Total number of samples*100), and species accumulation curves for all regions and moor types were created.

Sequencing summary
Overall read quality scores were lower for chip 3 compared to chips 1 and 2, with 56% of the reads being filtered out, compared to 17% and 38% for chips 1 and 2 respectively (Table S4). In total, after initial quality filtering by the PGM suite software, 2,037,175 raw reads were obtained from three Ion PGM chips. About 55% (1,126,242 reads) of these reads contained both a forward and a reverse barcode. Primer dimers could not be completely eliminated prior to sequencing due to the long length of the fusion primers (> 52 nts) and high concentrations of degenerate primers used, resulting in 37% of the raw reads with both barcodes being filtered out. Further, non-specific amplification (e.g. reads longer than the target) contributed to 26% of the raw reads with both barcodes. In total, 416,590 reads contained both a forward and reverse barcode and had the expected target size. Of these reads 78% (322,862 reads) were identified to a species level but the remaining 93,728 reads could not be identified with confidence to any taxonomic level using both the reference database and NCBI BLAST.
Two samples (Unmanaged moor England and Walkedup moor Wales) did not yield any prey reads with both COI and Cytb. Another two samples (Orkney Unmanaged moor and English DGM) only showed hen harrier (22 reads) with Cytb, and one sample (English DGM) contained 99.4% hen harrier reads and no prey OTUs with more than 0.5% relative abundance with Cytb (total read count in the sample -1401 reads). Mean read count per sample assigned to species with COI was 4433, with a median of 2067, ranging between 22 and 23903 reads. The number of OTUs per sample was not correlated with total read count per sample (Pearson correlation test, r = 0.17, p = 0.25, Fig. S2A). Mean reads per sample with Cytb was 1825, with a median of 915, ranging between 3 and 15,002 reads. The number of OTUs per sample detected with Cytb was also not correlated with the total read count per sample (r = 0.04, p = 0.78, Fig. S2B). Further, relative read abundance of hen harrier in the samples was not significantly correlated with the total read count per sample for both COI and Cytb (r = − 0.17, p = 0.24 and r = − 0.15, p = 0.29 respectively, Figs. S2C and S2D). In addition, there was no significant correlation between the relative read abundance of hen harrier in the sample and the number of OTUs identified in the sample (r = − 0.19, p = 0.19 for COI and, r = − 0.19, p = 0.19 for Cytb, Figs. S2E and S2F). Finally, relative read abundance of hen harrier showed a significantly positive correlation between both targets i.e. COI and Cytb (r = 0.44, p = 0.0012, Fig. S2G), and number of OTUs detected per sample was also found to have a significant positive correlation between both targets i.e. COI and cytb (r = 0.8, p = 3.8e−12, Fig. S2H).
A total of 310 OTUs (56 species) were identified with COI and 217 OTUs (43 species) were identified with Cytb (Fig. S2H). In both cases, the species accumulation curves did not reach the asymptote (Fig. 1). Combining COI and Cytb OTUs showed a higher detection of species (62 species) than both markers individually and also did not reach the asymptote. In samples where both targets yielded prey reads i.e. species were identified with both COI and Cytb, OTUs represented 37 species but COI more frequently detected these species than Cytb e.g. Skylark was identified in 30 samples with COI but Cytb showed presence only in one of these samples.

Control samples
PCR negatives amplified with COI did not yield reads identified to species, except 12 reads of Black grouse (Tetrao tetrix) on chip 2. PCR negatives amplified with Cytb contained exclusively human (Homo sapiens) sequences on chip 1 (6 reads) and chip 3 (48 reads), and no reads on chip 2. DNA extraction negative 1 (chip 1) contained sheep (Ovis aries) and two species of sturgeon (Acipenser ruthenus and Huso huso) (all used in the laboratory). Extraction negative 2 did not yield any reads on chip 2, however on chip 3, 22 reads of Meadow pipit were recovered with COI, but with Cytb, the same extraction negative resulted in human (129 reads), Rainbow trout (Oncorhynchus mykiss) (15 reads), and pig (Sus scrofa) (2 reads) (all used in the laboratory).
Positive controls showed that adding a blocking primer consistently increased reads assigned to potential prey species, especially in chips 2 and 3 where a more accurate Qubit quantification and a lower prey to predator ratio was used. Hen harrier was present up to 53% using COI and up to 26% using Cytb without blocking but never completely eliminated using Cytb with blocking (Table 1). Rabbit (Oryctolagus cuniculus) (COI) and Canada goose (Branta canadensis) (Cytb) had the least number of mismatches to the blocking primers, with only four differences, but results showed that these species showed no reduced amplification with introduction of blocking primers (Table 1). With COI, hen harrier was detected in 57 of 61 samples with mean reads of 16.1% (Range 0.02-96.2%) and a median of 1.95%. Predator blocking with the Cytb target resulted in hen harrier being detected in all but 1 sample with a mean of 38.3% (Range 0.26-100%) and a median of 32.33%.

Broods
From four broods, all the available chick samples were sequenced to see whether taking a single individual represented the total detectable prey diversity within the brood (Fig. 2). The three individuals from a brood on an English DGM revealed a total of five prey species. Individual 1 and 2 showed the same species, however, there were differences between the two markers in detection of meadow pipit and brown hare. Individual 3 of the same brood showed the addition of Common shrew (Sorex araneus), which was detected by both markers although Brown hare (Lepus europaeus) was not detected in Individual 3. The second brood from an English DGM shows a similar species composition between individuals, although, chicken (Gallus gallus) was detected in two of the four individuals, and Brown hare was detected in only one individual. Individual 4 of this brood yielded only hen harrier reads with Cytb. The brood from Orkney showed the highest number, with an overall 12 species. Individuals 1 and 4 showed the highest number, individuals 2 and 3 showed low overall read counts and presented only a small selection of prey species detected in the other two. The Scottish DGM brood showed three main species occurring  in all individuals and five species which were detected at low relative abundance only once in different individuals.

Overall prey species diversity
Of the 56 and 43 species identified with COI and Cytb respectively, 31 and 26 occurred only once or twice (Fig. 3). Both COI and Cytb combined showed 62 species from a wide range of taxonomic groups, predominantly avian and small mammalian species but also one amphibian, the common frog (Rana temporaria) detected in a single brood. All prey species that were recovered were native species occurring in Great Britain except Turkey (Meleagris gallopavo) which was present in one English DGM sample.
Other unusual domestic species detected were goat (Capra aegagrus) found in a single brood from an English DGM, and badger (Meles meles) in three broods from English DGMs across two different years (one brood in June 2006 and two broods in July 2007). Other raptor species were recovered, most notably red kite (Milvus milvus) in a brood from a walked-up moor in Wales, and in a brood from an unmanaged moor in mainland Scotland. Other raptors such as Buzzard (Buteo buteo), Kestrel (Falco tinnunculus), Merlin (Falco columbarius), and Sparrowhawk (Accipiter nisus) were also detected in one or two broods from walked-up moors in Wales. Besides raptors, two species of owls were also detected, Barn Owl (Tyto alba) was present in 6 broods (4 from walked up moors in Wales and 2 from Historical DGMs in mainland Scotland), and Tawny owl (Strix aluco) was detected in one brood from Wales. Finally, aquatic birds such as Red-throated diver (Gavia stellata) in a brood from Orkney, and Teal (Anas crecca) and Gadwall (Anas strepera) in a brood on an unmanaged mainland Scottish moor.

Regional and moor-type comparisons
Eight samples out of a total of 49 (16%) used in the analyses with COI and 11 samples out of a total of 50 (22%) used in the analyses with Cytb did not reach the 400 read threshold and were discarded for the moor-type comparison. For another five samples, precise moor-type was unknown and were therefore excluded from the moor-type comparison. Total number of broods used for the comparisons are shown in Table 2. We observed regional differences in prey With COI and Cytb, 5 and 4 species were found in more than 10 samples  Walked-up 7 6 6 6 6.7 5.3 8.2 6.5 5.0 7.5 5-9 3-9 6-12 + + + + − + species e.g. the endemic Orkney vole (Microtus arvalis orcadensis) was only detected in broods from Orkney and small birds of prey were only detected in Wales. The FOO (Fig. 4) shows that Meadow pipit, Red grouse, Wren (Troglodytes troglodytes) and Skylark (Alauda arvensis) occur in all regions and moor types. Field vole was present in all regions and moor types except Orkney and one brood from an unmanaged English moor. Brown rat (Rattus norvegicus) was only identified in samples from Orkney, and Rabbit, Starling (Sturnus vulgaris), and Pheasant (Phasianus colchicus) were abundant in Orkney but almost absent from other regions. None of the species detected were restricted to English DGMs except for four prey species that are not part of a hen harrier's natural diet in Great Britain: Chicken, Turkey, Badger and Goat. Northern wheatear (Oenanthe oenanthe) and Common shrew were the only prey species which were not detected on unmanaged moors but seen in DGMs and walked-up moors (Northern wheatear occurred once on an English DGM and once on a walked-up moor in England). Fourteen species were detected only in broods from walked-up moors and 24 species occurred only on unmanaged moorland. The FOO of Red grouse was 82% (9/11 broods) to 100% (8 broods) in English DGMs with COI and Cytb respectively. Red grouse FOO in unmanaged moors in mainland Scotland was 80% (4/5 broods) and 67% (2/3 broods) with COI and Cytb respectively but a lower FOO for Red grouse was seen in Orkney with 44% (4/9 broods) with both COI and Cytb, and in walked-up moors in Wales with 33% (2/6 broods) and 17% (1/6 broods) with COI and Cytb respectively (Fig. 4). The number of broods from mainland Scottish DGMs and Historical DGMs were too low for meaningful FOO estimates. A summary of the number of species per brood using both COI and Cytb based on region and moor-type is shown Fig. 4 The frequency of occurrence (FOO) separated based on region and moor type. The black bars are based on the detection with Cytb and the grey bars are with COI. The species order on the X-axis is based on the accumulative occurrence across the region and moor-type with the most commonly identified species in the first position on the left. Thick bars represent species that were identified with only one marker across the data set 1 3 in Table 2. Comparing moor types for which we had ≥ 6 broods, on average, the lowest number of species per brood was observed in samples from DGMs in England (Mean 3.5, median 3.5). Higher numbers of species per brood were recovered in unmanaged mainland Scottish moors (Mean 6, median 5) and unmanaged Orkney moors (mean 6, median 6). Broods from walked-up moors from Wales showed the highest number of species (mean 8.2, median 7.5). Although only one brood was available from both English walked up and unmanaged moors, 6 and 11 species respectively were detected, both higher than the mean seen in English and Scottish DGMs. Further, the brood from the unmanaged English moor showed 4 species more than the brood with the highest number of species recovered in English DGMs and 5 species more than the brood with the highest number of species recovered in Scottish DGMs ( Table 2). The two broods on a Historical DGM in mainland Scotland showed a high number of species (10 and 14), which is beyond the range seen in English DGMs but within or close to unmanaged moors in Orkney and mainland Scotland (Table 2).
Species accumulation curves (Fig. 5) show that in samples from DGMs in England a total of 13 species were recovered after 12 broods. Species richness from unmanaged moors in Orkney was at 21 after 12 broods. Further, there was no overlap in the 95% CI of English DGMs and unmanaged moors from Orkney (at sampling of 8 or more broods, upper 95% CI was 13.8 for English DGMs and lower 95% CI for Orkney was 14.4). The samples from walked-up moors in Wales show the highest species richness with 24 species after 6 broods. There was no overlap in the 95% CI for broods from walked-up moors from Wales and English DGMs at sampling of 3 or more broods (Lower 95% CI for Wales was 13.6 and higher 95% CI for English DGMs was 10.5). No overlap in the 95% CI was seen at 5 or more broods between broods from unmanaged moors in Orkney and walked-up moors in Wales (Lower 95% CI for Wales was 20.5 and higher 95% CI for Orkney unmanaged was 19.5). The 95% Confidence Interval (CI) for unmanaged moors in mainland Scotland was wide and overlapped with those from all other moor types. Samples from DGMs in mainland Scotland and historical DGMs in mainland Scotland resulted in 10 and 16 species respectively but numbers of samples from these moor types were low (n = 3 and 2 respectively).
A summary of the presence and absence of OTUs assigned to broad taxa is also shown in Table 2. Red Grouse, Passeriformes and Rodents/Eulipotyphla were present in all regions and moor types (only regions/ moor-types with > 4 broods considered). Lagomorphs were not detected in walked-up moors from Wales. Waders (Charadriformes) were detected in all regions except DGMs (Table 2).

Fig. 5
Species accumulation curves separated by moor type. Data from COI and Cytb were combined. Species presence was recorded when detected with at least one marker. Broods with < 400 reads were not used. Samples from English DGMs reached the asymptote after 12 samples with 13 species. Samples from unmanaged moors from Orkney showed 9 more species after the same number of samples as English DGMs. Samples from unmanaged moors from mainland Scotland showed the same species accumulation after 6 samples as Orkney after 6 samples although the 95% CI is wider. The historical DGMs showed double the number of species after a single sample as three from DGMs from mainland Scotland 1 3

Buccal swabs and dietary analyses
We have shown for the first time that non-invasive buccal swab samples from chicks can be used successfully for identification of prey species in birds. Buccal swab samples have several advantages over faeces, the most commonly used non-invasive sample for DNA based diet analysis e.g. (1) fewer PCR inhibitors are co-extracted with DNA, (2) there is less chance of contamination from environmental DNA, (3) DNA undergoes less degradation than in faecal samples which could have been deposited long before the sample was taken, and (4) the study species is directly known whereas identification of faecal source can sometimes be difficult (Morin et al. 2016). Success rates of recovering at least one prey species within a sample (97% with COI and 92% with Cytb) and those with > 400 reads used for region/moor-type comparisons (84% with COI and 78% with Cytb) in this study are similar to previous studies using faecal samples e.g. 95% in Louisiana Waterthrush (Parkesia motacilla) (Trevelline et al. 2016); 97% in Norway lemming (Lemmus lemmus) (Soininen et al. 2013). Other studies have reported far lower success rates from faeces e.g. 47% to 71% in Atlantic puffin (Fratercula arctica) (Bowser et al. 2013); 73% in Panthera species (Xiong et al. 2017).
Another bias found when using faecal samples is differential digestion e.g. Deagle and Tollit (2007) showed that proportion of prey reads per species was significantly different to expected values in captive Steller sea lions (Eumetopias jubatus). Minimal digestion occurs within the beak of a bird, so buccal swabs may be a more accurate representation of prey species for dietary studies over a longer time interval. With COI and Cytb, we identified on average 6.5 species per sample with a range between 1 and 15 species. Deagle et al. 2010 reported persistence of prey DNA in penguin faeces for at least 4 days after ingestion and it is possible that persistence is as long or longer in buccal swabs since 15 species were detected in 1 sample and chances of a chick being fed 15 species in 1 day is unlikely. Further, species accumulation curves (Fig. 1) suggested that not all the species were recovered with 51 buccal swabs i.e. asymptote was not reached. Although related to the range of prey species consumed by a predator i.e. asymptote likely to be reached sooner when range of prey species is low, these results highlight the usefulness of buccal swabs for dietary studies in avian species.
In-brood variation in prey detection was relatively high in some broods, with individuals in the same brood missing up to five species found in the complete brood (Fig. 2). Although this could be a result of the difference in read coverage between samples which was present in three of the four broods including samples with < 400 reads, this is unlikely because (i) the number of species per sample was found not to be correlated to the total read count per sample ( Fig. 1a and b), (ii) Orkney individual 2 had more than 400 reads for Cytb but showed lower number of species than individuals 1 and 4, and (iii) Scottish DGM individuals 2 and 3 had < 400 reads for Cytb and COI + Cytb respectively but showed presence of novel prey species not detected in individual 1. It is instead possible that differences in age of the chicks and their feeding behaviour resulted in this variation.

Primers, targets, controls, and predator blocking
Highly degenerate primers to COI and Cytb (384 and 648 oligos per primer) were needed to amplify across such a wide range of species. This resulted in a proportionally small percentage (29%) of the reads that could be identified to species due to presence of non-specific amplicons and residual primer dimers (Table S4). Using two target regions for a metabarcoding based diet study is not common. Other studies that have used two target regions also reported differences in recovery of prey species within the same samples with the two regions (Bowser et al. 2013;Gerwing et al. 2016). We selected COI because it is the 'Barcoding of Life' marker for animals (Hebert et al. 2003). Cytb was added as the second marker since it has been previously reported to resolve phylogenies better than COI in mammals (Tobe et al. 2010). Bowser et al. (2013) reported that that only 20% of all taxa were detected with both COI and 16S. In this study, difference between COI and Cytb was less apparent, with 59% of the taxa identified with both markers, 31% of the taxa only identified with COI and 10% of the taxa only identified with Cytb. Cytb identified on average 1 species less per sample than COI. Blocking of predator DNA was highly effective with COI but not as effective with Cytb possibly due to the amplification of a Nuclear Mitochondrial DNA segment (Numt) since a large proportion of the Cytb hen harrier reads matched the hen harrier Cytb sequence in the reference database and GenBank at only 92-93% homology (data not shown). No stop codons were found in these short sequences (data not shown) but multiple Numts have been previously reported from falcons (Nacer and do Amaral 2017). However, some species were more often only identified with Cytb e.g. pheasant was only identified using Cytb in seven samples and Orkney vole was identified in five samples using Cytb but twice with COI. This highlights the advantages of using two target regions for such studies, previously also mentioned by Bowser et al. (2013).
Although samples were pooled in equimolar concentrations, read coverage was variable within and between chips with both gel (chips 1 and 2) and Qubit quantification (chip 3). Such variability has also been reported in previous studies (Bowser et al. 2013). The addition of positive control samples on every chip gave useful insight into amplification of mixtures of DNA as well as into effectiveness of the blocking oligos. Results showed that all species were successfully recovered on all chips although differential levels of amplification were seen e.g. with COI, common shrew showed 9% of the reads where 20% was expected (Table 1). Adding blocking oligos resulted in an increase in detection of prey reads in all cases. The two species with the minimum number of mismatches to the blocking oligos, rabbit (COI) and Canada goose (Cytb) (Table S1), were both successfully amplified.
PCR negatives showed minimal contamination with prey species, with only 12 reads of Black grouse (Tetrao tetrix) which was abundant (546 reads) in one of the samples on the same chip. Extraction negatives did not contain prey species except 22 reads of Meadow pipit recovered in chip 3 (only 0.24% of the average reads per sample obtained on chip 3 and therefore < the 0.5% filtering applied to samples). The same extraction negative was used on chip 2 but did not yield any reads, suggesting that high coverage on chip 3 (318 chip) increased the chance of detecting contamination and low abundance OTUs. Other species in extraction negatives were all species handled in the lab such as sheep, sturgeon, trout, and pig. Since different sized chips were used, we used a 0.5% filtering criteria across all chips (Albaina et al. 2016 andAziz et al. 2017) to eliminate all low abundance OTUs. We also discarded samples with low overall read coverage (< 400) for region/moor type comparisons to ensure comparable data.

Overall prey species diversity
A large number of species (62) was identified in 51 broods using metabarcoding which is almost twice as high as with conventional studies (37 species) based on 2614 h of observation (Redpath and Thirgood 1997). Our results confirm that hen harrier diet is diverse. Previous studies have reported species from several different orders of birds, mammals, reptiles and amphibians (Balfour and Macdonald 1970;Redpath and Thirgood 1997). Former diet studies in the hen harrier (e.g. Redpath and Thirgood 1997;Redpath et al. 2001) have been limited to small geographical areas and no diet data from England and Wales were previously available. Using metabarcoding, our results show for the first time a snapshot of the diet of hen harriers across the whole of Great Britain. The majority of prey species occurred only once or twice across samples, with a few species being most frequent e.g. Meadow pipit, Red grouse, Wren, Skylark, and Voles (Field vole and Orkney vole). This is in agreement with previous diet studies from Scotland (Picozzi 1978;Redpath and Thirgood 1997). Birds of prey such as Kestrel (Falco tinnunculus), Merlin (Falco columbarius) and owls identified in this study, especially in Wales, have previously been identified by Balfour and Macdonald (1970) in Orkney. Hen harriers are known to be general and opportunistic hunters (Redpath and Thirgood 1997). This is supported by the results in this study, with a high number of species identified within a few broods.
More unusually we detected Red kite and buzzard which are much larger in size and may appear unlikely prey items. However, relative read abundance was 0.7, 4.2 and 10.8% for red kite in the three broods where this species was detected with COI in Wales and mainland Scotland and 14.1% for buzzard detected in walked-up moors from Wales. These percentages are unlikely to represent sampling of red kite/ buzzard saliva from a previously eaten carcass subsequently scavenged by a hen harrier but must instead represent rare cases of predation on red kite/buzzard chicks. Although Red kites nest in woodlands where hen harriers do not hunt, there are reports of Red kites nesting close to the edges of woods allowing a direct aerial route for these large winged adults (Carter and Grice 2000) and on occasion, perhaps offering predatory opportunity to a foraging hen harrier.

Regional and moor-type comparisons
All prey species detected occurred in the regions where broods were located e.g. Orkney vole was only detected in Orkney, Red-throated diver found in a brood from Orkney where large numbers of these birds are found, and brown rat being common and only identified in samples from Orkney. Also, species present in a large number of Orkney samples such as Rabbit, Starling, and Pheasant were rarely present in other regions (Fig. 4). These concur with previous diet studies e.g. Picozzi (1978) and Redpath and Thirgood (1997) did not record brown rat outside Orkney; Balfour and Macdonald (1970) and Picozzi (1980) identified up to 47% of the prey as Rabbit, 10% as Starling, 4% as Orkney vole in Orkney, but these were not commonly observed in mainland Scotland (Picozzi 1978;Redpath and Thirgood 1997). These previous studies had not identified pheasant as prey in Orkney but 6 out of 12 broods in this study were found to contain pheasant, most likely due to the more recent population expansion of pheasants in Orkney (Balmer et al. 2013) making them more available as prey species for hen harriers.
None of the natural prey species detected was restricted to DGMs, while 24 species occurred only on unmanaged moors from all regions and 14 species were only recovered from Walked-up moors in Wales. Chicken and Turkey were found in broods in English DGMs during a time when diversionary feeding was commonly undertaken. Whilst poultry chicks and white rats (Rattus norvegicus) are commonly used for diversionary feeding of hen harriers (Ludwig et al. 2018), turkey was detected once in this study at a relative 1 3 read abundance of 1.6% for COI, suggesting that on this occasion, turkey may have been used. One brood from an English DGM was also found to contain goat at 1.3% relative abundance for which there is no obvious explanation since goat meat is not commonly used for diversionary feeding. Badger was found over two consecutive years in three broods on English DGMs indicating that hen harriers possibly scavenged dead badger carcasses. Badgers are sometimes caught in snares used for predator control on DGMs (Thompson et al. 2016).
We found lower numbers of prey species identified per brood and a lower overall prey species richness in DGMs from England compared to unmanaged and walked-up moors (Figs. 4,5,and Table 2). Further, the species accumulation curve for English DGMs contains 4 species (badger, goat, turkey and chicken) that are not part of the natural diet, making the curve an overestimate of prey species richness. The alternative form of grouse shooting in the walked-up form practised in Wales, where hunters walk up moors with lower densities of grouse and flush the grouse with dogs ( Thompson et al. 2016), showed higher species richness, with large numbers of prey species detected. Broods from walked up moors in Wales showed a steeper increase in species accumulation than broods from unmanaged moors from Orkney due to the detection of many different types of raptors and owls in Wales which were not detected in Orkney.
The FOO of Red grouse was highest among broods from English DGMs but lower in broods from unmanaged moors in Orkney and walked-up moors in Wales (all with ≥ 6 broods for both COI and Cytb per region and moor-type, Fig. 4). Although the FOO of Red grouse was also high in broods from unmanaged moors in mainland Scotland (5 for COI and 3 for Cytb), we observed a high level of prey species richness in these broods, indicating that red grouse represents just one of many prey species available to hen harriers in these moors. In diversionary feeding experiments e.g. Ludwig et al. (2018), up to 76% of the diversionary food provided was taken by hen harriers, confirming that hen harriers are opportunistic and frequency-dependent predators, taking more of any species that is abundant or easily available. Although sample sizes per region and moor-type were too low for elaborate statistical analysis, a clear difference was visible in species richness and the number of unique species recovered. Both values were lower in English DGMs compared to unmanaged moors in Orkney and mainland Scotland as well as Walked up moors in Wales (all with ≥ 6 broods).
The two broods from a historical DGM in mainland Scotland revealed high numbers of prey species (14 and 10 respectively). At the time of sample collection, this DGM with low intensity management was surrounded by a mosaic of unmanaged moorland and other habitats such as young plantations and it is likely that these hen harriers hunted a broad variety of species in the mosaic of different habitats. Although DGMs have been reported to have two to five-fold higher breeding densities of Golden plover, Curlew and Lapwing than on unmanaged moorland (Tharme et al. 2001), we did not detect Golden plover in any brood, detected Lapwing once in Orkney where there are no DGMs, and detected Curlew in three broods, two from Orkney and one from the historical DGM in mainland Scotland.
In conclusion, our results show a high incidence of Red grouse and low prey species richness in the diet of hen harriers in DGMs. Additional studies are needed, but if our findings are due to an over-abundance of Red grouse in DGMs but no overall difference in prey species abundance and diversity across moor-types, a lowering of Red grouse density in DGMs could shift predation pressure from hen harriers off Red grouse to other available prey species. If, however, our findings indicate that the number of prey species available to hen harriers in DGMs is limited to Red grouse and a few other species because of lower prey species abundance and diversity due to intense management for Red grouse, the conservation conflict is likely to worsen in the future if management is further intensified.