Introduction

Insect-eating birds across the globe consume 400–500 million tons of insects annually (Nyffeler et al., 2018), surpassing the 350 million tons of meat that humans consume each year (Hicks et al., 2018). However, insectivorous bird populations are particularly at risk due to climate change. In North America, some terrestrial insectivorous bird populations have decreased by 33% in the past five decades (Tallamy & Shriver, 2021). Similarly, many insectivorous and non-insectivorous bird populations have faced a substantial decline over the past 30 years in Europe (Inger et al., 2015).

Upper elevation birds are particularly at risk because mountains are expected to be more affected by climate change than lowland areas due to faster and enhanced warming (Mountain Research Initiative EDW Working Group, 2015). For example, high-elevation populations of Canada jays (Perisoreus canadensis (L., 1766)) declined 50% over a 30-year period, and the decline was attributed to warmer and more variable weather (Sutton et al., 2021). This change in weather pattern increased the number of freeze–thaw events, which caused an increase in the spoilage of cached food items. The survival of another high-elevation Passerine, the white winged snowfinch (Montifringilla nivalis (L., 1766)) is in doubt because its foraging behavior is closely tied to snow retreat conditions which are becoming increasingly less consistent (Resano-Mayor et al., 2019). Finally, Barras et al. (2021) found that elevated ambient temperatures at the tree-line in the Swiss alps negatively affected nestling prey provisioning rates of the Alpine ring ouzel (Turdus torquatus alpestris L., 1758).

The flora and fauna of the Pyrenees mountains are especially threatened due to climate change and land use shifts (OPCC-CTP, 2018). The snowpack is warmer than many other mountain ranges across the world and thus particularly sensitive to slight changes in ambient temperature (Lopez-Moreno et al., 2017), and the decline of agropastoral practices in the Pyrenees has led to transitions of open grassland into forest (Roura et al., 2005). The Pyrenean tree-line is also shifting upward (Ameztegui et al., 2016). Tree-line dynamics globally are affected by a variety of factors, including precipitation, tree community composition, and soil structure (Grace et al., 2002; Körner, 2012). The upward shift and densification of the tree-line in the Pyrenees is generally linked to local agricultural abandonment (Batllori & Gutiérrez, 2008), but there can be locally important factors such as slope morphometry and lithology (Feuillet et al., 2020).

It is within this context that we examined the diet of Passerines at elevations located below and at tree-line. The diet of many European Passerine birds, e.g., Paridae, has been examined closely, even though most studies were limited to estimating the diet of nestlings using either methods that are invasive (neck collars (Barba & Gil-Delgado, 1990; Pagani-Núñez et al., 2011) and stomach flushing (Senécal et al., 2021)), noninvasive but less detailed (e.g., cameras and stable isotope analysis (Currie et al., 1996; Şekercioğlu et al., 2023, respectively)), or lethal (gizzard extraction (Sehhatisabet et al., 2008)). Recent advances in DNA metabarcoding technology (hereafter metabarcoding) have increased our ability to analyze the diet of adult Passerines in a noninvasive manner at a high level of taxonomic classification (see Crisol-Martínez et al., 2016; Ribeiro et al., 2019; da Silva et al., 2020; Shutt et al., 2020).

Metabarcoding supports high-throughput (i.e., massively parallel) taxonomic classification within a sample (Bush et al., 2019). A short portion of a gene (barcode) from an environmental or biological sample is amplified by a primer designed to provide taxonomic resolution of a target organism or taxonomic group (Deagle et al., 2014; Hajibabaei et al., 2007). However, there are many fundamental limitations of fecal metabarcoding. For example, raw prey abundances cannot be determined from the number of reads in a similar DNA sequence, and relative prey abundance is difficult to recover because of technological and biological biases including primer mismatch and differences in PCR amplification due to primer sequence length (Deagle et al., 2013; Krehenwinkel et al., 2017; Piñol et al., 2015). Sample contamination and differing rates of DNA preservation in the gut can also present issues (Galan et al., 2018; Nielsen et al., 2018).

A handful of studies have examined the diet of multiple adult Passerines species using metabarcoding (see Crisol-Martínez et al., 2016, Sottas et al., 2020, and Garfinkel et al., 2022). Most Passerine dietary metabarcoding studies have examined the adult diet of one species (McClenaghan et al., 2019; Moran et al., 2019; Shutt et al., 2020, and Snider et al., 2021 among others) or that of nestlings (see Rytkönen et al., 2019, da Silva et al., 2020, and Jarrett et al., 2020). The goal of our study was to examine the diets of co-occurring high-elevation Pyrenean Passerines and to assess the effects of open space, elevation, and season on them. We expected higher niche differentiation (i.e., difference in diet composition) in morphologically and behaviorally similar species and higher dietary richness and higher dietary overlap among species as spring progressed to autumn. Passerines that have similar traits often have competition-driven niche separation (Alatalo et al., 1986; Cowie & Hinsley, 1988; Sottas et al., 2020), and higher abundance of prey is linked to less dietary partitioning (Davies et al., 2022). In some species, we expected diet richness to positively correlate with the percentage of open space because patchier habitats have been shown to benefit some species but not others (Suarez-Seoane et al., 2002). Finally, we expected higher overlap in below tree-line plots because we expected the conditions to be more favorable to Passerines. Higher overlap is common in more favorable habitats (Hou et al., 2021).

Methods

Study area and feces collection

Ten plots were selected in grid format within a black pine forest (Pinus mugo species complex, hereafter P. mugo) in Vall d’Ordino, a valley located within three km of Vall de Sorteny Natural Park in the parish of Ordino, Andorra. Plots were situated between 1729 and 2352 m above sea level (hereafter MASL). Percent open space surrounding each plot (1000 m radius) was calculated using QGIS3.4 and the MCSA 2012 landcover map downloaded from the Institute of Andorran Studies (Centre de Biodiversitat de l'Institut d'Estudis Andorrans, 2012). Testing various scales of open space was out of the scope of this study; therefore, a landscape level scale of 1000 m radius was chosen as a 1000 m radius is a generally accepted radius for a landscape level measure of “forest amount” (Shoffner et al., 2018). Plots were characterized as “below tree-line” or “at tree-line” depending on positioning above or below the median elevation of all plots (i.e., 2077 MASL). In most of Andorra, tree-line occurs between 2200 and 2400 MASL and in a few areas between 2100 and 2500 MASL (Carreras et al., 1996).

Birds were captured using Ecotone mist-nets (9 m and 6 m long and 2.5 m high, with 5 shelves and a mesh size of 16 mm2) stretched between 4 m poles inserted perpendicularly in the ground. Three nets were used in eight of the sites, and two nets were used in the remaining two sites due to difficultly of access, and openness of site (see Supplementary material 1). Mist nets were deployed between May 15, 2018, and September 30, 2018, and number of net days for each site are also documented in Supplementary material 1. The start date for the field component of this study coincided with the date when snow historically has retreated from the Andorran tree-line. Mist nets were not set at a plot when there was precipitation or high winds to ensure good capture conditions and welfare of the bird. Once a bird was captured, it was placed in a single use individual paper bag. After defecation, feces were carefully removed from the paper bag using a single-use toothpick and stored in plastic vials. Vials were placed on ice in the field and transferred to long-term deep-freeze as soon as possible. All birds caught were identified at species level, ringed, aged, sexed, and measured following standard ringing procedure. Birds were handled by certified ringers and all the procedures approved by the Environment and Sustainability Department of Andorran Government.

A total of 132 fecal samples were collected. No bird was captured twice. Samples collected in May and June were considered to be from the spring, July and August samples were considered to be from the summer, and September samples were considered to be from the autumn.

DNA extraction and amplification

DNA in the collected fecal samples and in four DNA extraction blanks (i.e., vials with no fecal samples that served as contaminant control) was extracted using the QIAamp DNA Stool Mini Kit (QIAGEN) following the manufacture’s protocol with modifications suggested by Davies et al. (2022). Only samples of similar size (approximately 3 mg) were processed. Three samples were smaller than 3 mg and thus discarded from analysis. DNA concentration was quantified using a Qubit 3.0 Fluorometer (Thermo Fisher Scientific Inc.). Amplicon library preparation and sequencing were carried out by the Georgia Genomics and Bioinformatics Core (University of Georgia, Athens GA, USA) as follows: the concentration and size distribution of the DNA were determined using SYBR fluorometry and a fragment analyzer, respectively. The DNA samples were normalized to the same concentration before library preparation. The amplicon library was prepared beginning with a 2X KAPA HiFi HotStart ReadyMix (Roche Sequencing & Life Sciences). The first PCR was assembled with target-specific primers. Initial plans called for the use of a primer amplifying a longer region, but a preliminary experiment (data not shown) indicated higher efficacy of a shorter primer, the mini-barcode mitochondrial primer (ANML). ANML amplifies a smaller 180 bp segment on the cytochrome oxidase C subunit 1 (COI) (Jusino et al., 2017). In the second round of PCR, unique Illumina Indexed primers (i5 and i7) were added to each sample and at the end of each round of PCR, amplicon was purified with bead cleanup. Final libraries underwent quality checks by SYBR fluorometry, Fragment Analyzer and qPCR (KAPA Library Quantification—Roche Sequencing & Life Sciences). Libraries were normalized and pooled equally. The pooled amplicons were sequenced on the MiSeq platform (PE250) (Illumina).

Bioinformatic analysis

Within the QIIME 2 2020.6 environment, tagged feces sequence reads generated from the Illumina MiSeq sequencer were demultiplexed and primers trimmed to create fastq files (Bolyen et al., 2019). We then used the DADA2 pipeline for further downstream analysis, which created a table of amplicon sequence variants (ASVs) rather than traditional operational taxonomic units (OTUs), thereby improving reproducibility, comprehensiveness, and accuracy (Callahan et al., 2016). Potential contaminants in the ASV table using the “frequency” method were identified by the package Decontam using the default threshold of 0.1 (Davis et al., 2018). The “frequency” method identifies contaminants whose concentrations in blank samples are inversely correlated with sample DNA concentrations, based on a user-specified threshold. As the metabarcoding workflow introduces quantitative bias into results, ASV raw counts were transformed into a presence/absence matrix (Martoni et al., 2022).

ASVs were taxonomically classified by aligning sequences to those within the arthropod training database “tidybug” (O’Rourke et al., 2020) via the classy-sklearn naïve Bayes method implemented in QIIME 2’s q2-feature-classifier plugin. The full QIIME script can be found in Supplementary material 2. Each taxonomic assignment was further examined individually using the following protocol: (1) The geographic range of the assignment was assessed, and species that do not occur in Europe were removed. (2) Species considered to be rare in Europe but are not known to occur in the Pyrenees were flagged. (3) Flagged assignments were verified by submitting query sequences to the NCBI BLAST tool and assignments that did not score at or above 98% identity were removed. All single flagged species composed of multiple ASV sequences were aligned to check for sequencing error, and sequences above an 80% sequencing error were kept in the analysis (Brandt et al., 2021; Ritter et al., 2022). (4) If an ASV showed multiple hits with the same max score on the NCBI BLAST tool, the ASV was removed from the analysis.

Statistical analyses

A Wilcoxon–Mann–Whitney test was performed to test for correlations between open space and elevation as a discrete variable, using the percent of open space as a dependent variable. We calculated mean prey richness of all bird species and of the top two commonly collected bird species (Periparus ater (L., 1758) and Lophophanes cristatus (L., 1758)). Predictive roles of independent factors affecting of the combined prey richness of these two bird species (season, elevation as a continuous gradient, and open space) were found by fitting data to a negative binomial or Poisson model (GLMM) using the lme4 v.26 package (Bates et al., 2015) in R Version 1.3.1056 (R Core Team, 2021). Models were chosen through a combination of (1) residual plotting with the DHARMa package (Hartig, 2022); (2) model performance testing using Pearson, Kendall, and Spearman correlation coefficients; and (3) model accuracy evaluation by measuring the root mean square error and the mean-absolute deviation of each model. Plot was used as a random variable. We also tested “individual bird” as the only random factor and it yielded very similar results, quantitatively and qualitatively. Nevertheless, the “individual bird” choice also caused some models to fail convergence and these problems were more severe when both “plot” and “individual” were chosen together as a random factor. Therefore, our final choice for random factor was “plot.” Some models, however, did not accept any random factor and were thus run without one after numerical testing showed little differences between models with random factors and those without. After the prey richness model of each bird group was fit, post hoc Tukey tests were carried out to investigate error rates of the categorical factor of season. Comparison of beta diversity (i.e., dissimilarity of diet) among and within bird species was determined by a Jaccard dissimilarity matrix using the R vegan package, with a value of 1 indicating there were no shared species and a value of 0 indicating complete sharing of species (Oksanen et al., 2020). A Jaccard dissimilarity matrix was also used to examine the dissimilarity of diet of all individual birds combined and the individuals of the two most common bird species combined in regard to season and the discrete variable of elevation. Frequency of occurrence (number of times a prey item appeared in a fecal sample, divided by total number of samples) was calculated for the six most frequently occurring prey groups (Arachnida, Diptera, Coleoptera, Hemiptera, Hymenoptera, and Lepidoptera). Predictive roles affecting species richness within each prey group (season, elevation as a continuous gradient, and open space) were calculated using the same model selection process as bird groups. Sample coverage was examined using the iNEXT package to create sample size-based rarefaction curves and extrapolation curves (Chao & Jost, 2012; Chao et al., 2014).

Results

A total of 8.95 million sequence reads were produced in the feces samples, with ASV counts per feces sample ranging from four to 268,999 (Supplementary material 3). ASVs were taxonomically classified as the MOTUs of 713 arthropod classifications (representing six classes, 30 orders, 163 families, 466 genera) and then only MOTUs classified to genus or species were kept for a total of 590 arthropod MOTUs. We refer to these MOTUs as species.

Inter- and intra-species dietary richness of Passerines

Fecal samples were collected and analyzed from 14 bird species (Supplementary material 4). The bird species P. ater (coal tit), L. cristatus (crested tit), and Prunella modularis (L., 1758) (dunnock) accounted for 40%, 15%, and 14% (respectively) of the bird species from which samples were collected. No significant seasonal differences were found when all bird species were combined or when just P. ater and L. cristatus were combined (P values in Supplementary material 5). No captures of P. modularis were made in autumn. Seventy-three and 59 birds were caught in the high-elevation and low-elevation plots, respectively.

The mean prey richness per bird capture was 13.5 ± 6.9 species. Phoenicurus ochruros (Gmelin, 1774) (black redstart) displayed the highest species richness per bird capture, followed by Sylvia atricapilla (L., 1758) (blackcap) and P. modularis, although standard deviation bars do not show differences within these three species (Fig. 1).

Fig. 1
figure 1

Prey species richness per bird capture with standard deviation bars. Values in parentheses represent the captures (i.e., number of feces collected) of each bird species

Full size image

According to Wilcoxon–Mann–Whitney tests, the percent of open space was significantly higher in the plots at tree-line (mean: 88 ± 13.7) than the plots below tree-line (mean: 45.2 ± 17.5) (Z = 4.25, P < 0.001). However, GLMM results showed that the percentage of open space was not correlated with prey richness in the diets of P. ater and L. cristatus when combined, or when all bird species were combined (Beta estimates and P values in Supplementary material 5). There were no elevational differences in prey richness when all bird species were combined, or when P. ater and L. cristatus were combined (Supplementary material 5).

Inter- and intra-species dietary overlap of Passerines

Jaccard dissimilarity index showed very little overlap in the diet between and within bird species, and there was a mean dissimilarity of 0.91 ± 0.05 in the diet of the 14 bird species. Beta diversity was very high among species. Compared to each other, P. modularis /L. cristatus and P. modularis /P. ater had a higher rate of dissimilarity (0.86 and 0.83, respectively) than the rate between P. ater/ L. cristatus (0.73). The dietary variability within P. ater, L. cristatus, and P. modularis was high as well: P. ater (0.82 ± 0.02), L. cristatus (0.93 ± 0.05), and P. modularis (0.95 ± 0.04). Accordingly, sample size-based rarefaction curves indicated it would be necessary to capture over 100 more P. ater individuals than P. modularis and L. cristatus to reach 99% sample coverage (Fig. 2a). At 99% coverage, P. modularis is expected to have a higher diversity of diet than P. ater and P. cristatus, while P. ater is expected to have the lowest (Fig. 2b). The mean overlap among all individual birds captured was not different between the three seasons (spring, 0.9 ± 0.05; summer, 0.93 ± 0.07; autumn, 0.93 ± 0.06) (Supplementary material 6a). Similarly, when the mean overlap between the two most captured birds (P. ater and L. cristatus) was calculated by season, no difference in overlap was recorded (spring, 0.92 ± 0.07; summer, 0.91 ± 0.08; autumn, 0.90 ± 0.07) (Supplementary material 6a). The mean overlap between all birds and between P. ater and P. cristatus in the plots at tree-line and below tree-line was similarly low (see Supplementary material 6b).

Fig. 2
figure 2

Sample size-based rarefaction curves (solid line) and extrapolation curves (dotted line) with 95% confidence intervals. a Number of birds caught per sample coverage. Numbers in parentheses indicate bird species and number of individual birds caught at 99% sample coverage. To reach 99% coverage, it would be necessary to capture over 100 more P. ater individuals than P. modularis and L. cristatus. b Prey diversity of bird species per number of individual birds sampled. Numbers in parentheses indicate number of captures of birds per species necessary to reach 99% coverage, and prey diversity at 99% coverage with 95% confidence intervals. Prunella modularis is expected to have a higher diversity of diet at 99% sample coverage, while P. ater is expected to have lowest diversity

Full size image

Presence of prey species and prey species trends

Most prey species were rare; 58.64% of the prey species were collected only once (i.e., collected in one sample) (Table 1). However, eleven species were present in over 15% of samples (Supplementary material 7). Of these eleven species, five were conifer pests. Aphid conifer pests (Hemiptera) were the two species most likely to be present (Supplementary material 7). Diptera and Lepidoptera represented 21.78% and 20.11%, respectively, of all prey species (Fig. 3). GLMM results indicated that the richness of some groups of prey (Coleoptera, Hemiptera, and Hymenoptera) significantly varied between some seasons, and elevation and open space did not drive richness of prey groups (Beta estimates, P values, and seasonality Tukey tests in Supplementary material 8).

Table 1 Prey species distribution among bird feces collections
Full size table
Fig. 3
figure 3

Pie graph of relative species richness by order, i.e., percentage of total number of prey species belonging to each order. Only orders with relative species richness above 0.00% are included in pie graph

Full size image

Discussion

The difficulty in accurately and directly identifying species-level dietary components of adult insectivorous Passerines in a non-lethal manner is a quandary that has long confounded ornithologists. A metabarcoding approach allows for the direct study of species composition of Passerine feces. With this technique, we were able to determine that: (1) there was extremely high inter- and intra-species variability and low inter- and intra-species overlap in the diet of captured Passerines, a result that contrasts with some traditional Passerine nestling studies and agrees with many metabarcoding studies; (2) dietary richness did not correlate with season nor with open space, and there was no difference in dietary overlap relative to open space or elevation, which suggests that high-elevation Passerine birds in our study have high dietary mobility; and (3) Passerine diets were dominated by conifer pests.

There was high biological richness in our analysis of the 132 fecal samples: over 594 taxonomically classified arthropod species were identified. Most prey species were rare; 58.64% of the prey species (or 346) were recorded in only one feces collection. A recent metabarcoding study examining diets of 25 North American grassland Passerines also found very diverse diets and rare prey items among bird species and linked this inter-species dietary diversity to opportunistic availability (Garfinkel et al., 2022). Similarly, authors of other metabarcoding studies examining (respectively) Eurasian blue tits (Cyanistes caeruleus (Shutt et al., 2020)), Rufous hummingbirds (Selasphorus rufus (Gmelin 1788)) (Moran et al., 2019), and barn swallows (Hirundo rustica (L. 1758)) (McClenaghan et al., 2019) postulated that the high intraspecific dietary variation found in the birds’ diets was likely due to prey availability and dietary flexibility.

Despite a mean richness per bird capture of only 13.5 ± 6.9 species, the dietary overlap among and within bird species was very low. We expected lower overlap (higher niche differentiation) between closely related bird species, as it is well established that segregated foraging behavior occurs between closely related European insectivorous Passerines. When a potential niche is left unoccupied by a Passerine bird species, the species that most resembles the absent species in body size is the species most likely to fill it (for a review of geographic niche changes in insectivorous Passerines, see Alatalo et al. (1986)). Segregated foraging behavior makes sense in light of a study showing negative impacts upon a less dominant but closely related species sharing geographic space; Parus major (L. 1758) (great tit) nestlings raised sympatrically with C. caeruleus weighed less than those raised allopatrically, suggesting that a large overlap of resource utilization exists between the two closely related species (Torok & Tóth, 1999). Most dietary studies comparing insectivorous hole-nesting Passerines (mainly some combination of P. ater, P. major, L. cristatus, and C. caeruleus) have historically examined nestling diets and reported a high overlap when dietary components are classified to a combination of class and family (Grzędzicka, 2018; Michalski et al., 2011; Nour et al., 1998). One study that classified Passerine prey of Lepidoptera and Arachnida to species also found high overlap (Atiénzar et al., 2013). However, at least two studies have shown that P. major and C. caeruleus feed differing sizes of caterpillar prey to nestlings (Ceia et al., 2016; Torok & Tóth, 1999), a result that would not be observable in a DNA-based study such as ours. In Ceia et al. (2016) the authors showed that the composition of prey (classified to family) between the two bird species were similar, and they postulated that difference in prey size resulted from either interspecific competition between P. major and C. caeruleus or the segregation of bird foraging guilds; Cyanistes caeruleus is primarily a foliage-gleaner and more likely to come in contact with smaller instar caterpillars, while P. major are bark–foliage gleaners and therefore more likely to come in contact with later instar caterpillars. Numerous dietary metabarcoding studies, however, have found lower inter- and intra-species overlap, including a dietary metabarcoding study of insectivorous Malaysian Babbler species (genera within the families Pellorneidae and Timaliidae) which found lower than 35% overlap in three of the ten species (Mansor et al., 2022). Lower inter- and intra-species overlap was also found in two other metabarcoding studies that examined dietary overlap within European insectivorous Passerines species (see Rytkönen et al., 2019; Shutt et al., 2020).

While inter- and intra-species overlap was very low in our study, the Jaccard dissimilarity indexes displayed slightly higher index values between P. modularis/L. cristatus and P. modularis/P. ater than between the more closely related P. ater and L. cristatus. These results are likely due to diverging foraging habits, a finding supported by a metabarcoding study that found differing dietary compositions of ground- vs arboreal-foraging bird species in macadamia orchards (Crisol-Martínez et al., 2016). Periparus ater and L. cristatus both forage in trees (Alatalo, 1981; Hartley, 1987; Lens, 1996), while P. modularis are mainly ground feeders (Bishton, 1986).

Historical data report divergent timing of clutch-laying of closely related insectivorous hole-nesting Passerines (Sanz et al., 2010), and a more recent study reported that many Passerine nestlings are provisioned with differing prey types depending on the nestling’s development stage (Orłowski et al., 2017). Historical data also suggest that resident insectivorous Passerines in many temperate forests are less segregated in both foraging sites and dietary components in summer, when insect prey is more abundant. Insect prey in autumn and winter is less abundant, leading to resource partitioning and inter-species competition (Betts, 1955; Gibb, 1954; Lister, 1980). There has been disagreement over seasonal segregation and diet. For example, Ulfstrand (1977) found spatial segregation in summer compared to autumn, whereas Wagner (1981) and Almeida and Granadeiro (2000) found no significant seasonal spatial differences. Obeso (1987) found no spatial difference but did find significantly different dietary components. Finally, a recent study found high dietary overlap between communities of insectivorous Passerines during times of limited insect availability, in contrast to many studies that show high dietary overlap during times of high resource availability. The researchers postulated this result indicated that during times of very low food availability, bird species were unable to avoid competition (Kent et al., 2022). Our study found no seasonal differences in dietary overlap between the two most commonly captured birds (P. ater and L. cristatus).

Our results also did not show seasonal dietary richness differences of these two bird species combined or when dietary richness of all individual birds captured were combined (Supplementary material 5). Our seasonal richness results are different from the findings of an adult metabarcoding study done in a deciduous forest (Shutt et al., 2020). In that study, Shutt et al. (2020) tested 793 fecal samples collected from March through early May of C. caeruleus and linked rising seasonal dietary richness to rising herbivorous insect abundance and availability. It is possible that the lack of differences in seasonal dietary richness found in our study is a result of our sampling time frame (mid-May through late September) and/or the relatively lower number of samples taken in our study (132). Regardless, our results (lack of seasonal prey overlap among bird species, no difference in dietary richness from spring to autumn, and no difference in richness levels among all bird species) likely indicate that the birds captured in our study have a high level of dietary flexibility.

We expected dietary richness levels of ground and shrub foraging species, such as P. modularis, to be correlated with open space, as structure and composition of vegetation can be very influential in nestling success of some European Passerine species (Orłowski et al., 2017). We also expected higher dietary overlap in plots below tree-line; while response of insects to elevation is species specific (Hodkinson, 2005), many montane fauna either decrease with elevation or have a humped shaped distribution along an elevational gradient (McCain, 2009; Rahbek, 2005). Temperature is known to be a major driver in insect community structure (Bale et al., 2002), and temperature swings are wider at higher elevations in the Pyrenees (Navarro-Serrano et al., 2020). Therefore, insects are likely more abundant in plots below tree-line, and thus, insectivorous Passerines might be less segregated at these lower elevations. However, in our study percent of open space and/or elevation had no effect upon the richness of the bird diets, or when the most common birds were examined separately. While more birds were caught in plots below tree-line than in plots at tree-line, the diet composition within both these groups showed low overlap, i.e., the diet among birds in plots at tree-line had as much overlap as the diet of birds in plots below tree-line. More data would be needed to document and compare the diet of each 14 bird species we studied, but this lack of link between open space and elevation may indicate high mobility of the more common generalists that occupy high-elevation Pyrenean landscapes.

Our adult Passerine barcoding study was performed in P. mugo forests, so it is unsurprising that five of the eleven species most likely to be present in the bird feces from this habitat were conifer pests. The two species most likely to be present in the bird feces were conifer aphids; 53.03% of the samples contained Cinara pini (L., 1758) and 37.88% contained Eulachnus rileyi (Williams, 1911) (Supplementary material 6). Cinara pini is a common and native conifer pest in Europe. Eulachnus riley, however, is considered rare to uncommon in its native range in Europe and is considered a pest outside its native range (Blackman & Eastop, 1994), so it is interesting to find this species to be common in our study. Even though, as a group, the European Paridae are some of most intensely studied birds in the world (Betts, 1955; Cowie & Hinsley, 1988; Gibb, 1954; Lack, 1964; Ulfstrand, 1976), data are limited relative to their adult diets. However, the abundance of aphids in adult diets in our study is similar to two other studies: Shutt et al. (2020) reported aphids comprised three of the top ten prey taxa and had the highest presence incidences, and Betts (1955) found aphids comprised over 50% of the adult diet of three species of hole-nesting European Paridae in June. Birds play important roles in top-down control of forest arthropod populations (Fayt et al., 2005; Gunnarsson, 1995, 1996; Philpott et al., 2004; Schwenk et al., 2010). It is possible the Passerines in our study are shaping arthropod communities and causing a trophic cascade by affecting tree growth. Research examining trophic cascade affects by bird predation, however, have revealed complex interactions or mixed results (Gruner, 2004; Schwenk et al., 2010).

Regardless, we now know many high-elevation species are under pressure (Öztürk et al., 2015), and alpine-breeding Passerines such as Anthus spinoletta (L., 1758) (water pipit) are precipitously declining in some areas (Flousek et al., 2015). We caught only one specimen of A. spinoletta; the remainder of our species are elevational generalists and not confined to high elevation. While our data does not shed light on alpine specialists, the lack of differences in dietary overlap and diversity relative to open space and elevation, not to mention the extremely low levels of intra-species dietary overlap, suggest that adult diet may not be a constraining factor in populational growth of some generalist insectivorous Passerines in and around the Pyrenean tree-line. At least one European generalist insectivorous Passerine bird seems to display extreme plasticity in timing of egg-laying (Wesołowski et al., 2016), and other ecological requirements such as suitable nesting sites and the provisioning needs for nestling could be more plausible population constraints. The upward migration of the tree-line in the Pyrenees, a phenomenon likely caused by both land-use shifts and climate change (Batllori & Gutiérrez, 2008; Batllori et al., 2010), may therefore not be a significant factor affecting the diet of some adult generalist insectivorous Passerines.

Much remains to be discovered regarding the diets of adult European insectivorous Passerines (Cholewa & Wesołowski, 2011). In the future, metabarcoding will undoubtably continue to elucidate the relationship between birds, insects, and landscape and has the potential to reveal vast quantities of dietary data. Our results showed very high prey diversity and very little overlap within and among insectivorous Passerines. Spatial trends (open space and elevation) had little effect on prey diversity and overlap. While these data indicate that the dietary plasticity of the more common birds is high, more studies are needed to reveal dietary components of mountain species, such as A. spinoletta.