Multi-dimensional niche differentiation of two sympatric breeding secondary cave-nesting birds in Northeast China using DNA metabarcoding

Niche theory predicts that ecologically similar sympatric species should show differentiation in at least one of the main niche dimensions (time, space, and/or food). Here, we combined observations of breeding timing, nest site selection, and diet (the latter determined using DNA metabarcoding) to analyze the niche overlap and differentiation between two sympatric secondary cavity-nesting birds, the Japanese Tit Parus minor and the Yellow-rumped Flycatcher Ficedula zanthopygia. The results showed that 1) there were significant differences in the first egg laying date, length of the egg laying period, incubation date, and hatching date between tits and flycatchers, and the breeding time of flycatchers peaked later (about 30 days) than that of tits; 2) the two species had a large overlap in nest site selection, although the canopy coverage and shrub density of flycatchers were significantly higher than those of tits; and 3) the niche overlap in diet was minimal, with both species heavily relying on Lepidoptera (39.6% and 63.7% for tits and flycatchers, respectively), but with flycatchers consuming significantly higher percentages of Lepidoptera, Diptera, and Coleoptera than tits. The results indicate that these two sympatric secondary cavity-nesting species have significant niche differentiation in breeding time and diet, but little differentiation in nest site selection.

Methods
Study area

The study was conducted during two breeding seasons (2021–2022) in Zuojia Nature Reserve, Jilin, northeastern China (44°1ʹ–45°0ʹN, 126°0ʹ–126°8ʹE). The reserve has an area of 60 km². The altitude is 200–554.6 m, with a continental monsoon climate in the temperate zone. The average annual precipitation is 680 mm, and the average temperature is 5–6°C. The vegetation is a natural secondary broad-leaved forest (Deng et al., 2011), including Betula dahurica, Fraxinus mandschurica, Quercus mongolicus, and Tilia mandshurica. Artificial nesting boxes (12 cm × 12 cm × 25 cm inside, 4.5 cm diameter of the hole) were hung in the study area to attract the tits and flycatchers. The nesting boxes were randomly attached to trees at a height of about 2.5–3 m above the ground, with a minimum distance of 30 m between each box. The tree species and nesting orientation were assigned at random, and the number of nesting boxes was maintained at about 450 year round (Yu et al., 2017).

Data collection

Collection of breeding time data

We visited the nest boxes every 3 days after mid-March to monitor the reproductive activities and clutch sizes of these two species (Fan et al., 2021). We recorded the first egg-laying date, incubation date, and hatching date of the tits and flycatchers. The breeding time did not include the nestling period in this study due to the difficulty of accurately monitoring the fledging date of each nest. Breeding time data were continuously recorded for 152 and 26 nests of tits and flycatchers, respectively. In addition, to measure reproductive fitness, we recorded the brood size of both species and calculated the hatching rate (i.e., the ratio of brood size to clutch size). We collected hatching data from 57 nests of tits and 26 nests of flycatchers between mid-May and late June, when both species were breeding and experiencing interspecific competition.

Collection of nest site data

We measured 13 nest site characteristics (Li et al., 2023) of the nest boxes occupied by tits and flycatchers. The nest height (m) (NH) above ground, the nest tree diameter (cm) at breast height (DBH), and the average diameter (cm) at breast height (ADBH) of 10 trees within a 10 m radius circle centered on the nesting tree were measured using a tape measure. The nest tree height (m) (TH), average height (m) of 10 trees (ATH) within a 10 m radius circle centered on the nesting tree, the average height (cm) of 10 shrubs (ASH), and the shrub density (SD) within a 1 m radius circle centered on the nesting tree were measured using a laser rangefinder. We recorded the nest tree species (TS), the number of tree species (NTS), and the number of trees (NT) within a 10 m radius circle centered on the nesting tree and assessed canopy cover (%) (CC) above the nest tree. The CC was measured by standing 1 m away from the nest tree in each of the four cardinal directions and observing it vertically to assess it (Li et al., 2023). The entrance inclination (EI) of nest boxes was measured by a digital protractor (Weidu Electronics Co., LTD, Wenzhou, China). The orientation of the nest box entrance (OE) was measured by a compass. We collected data for 300 nests of tits and 34 nests of flycatchers. The nest site data included nests with complete breeding time data, as well as those with incomplete data (i.e., missing hatching dates or nests preyed upon by predators).

Collection of adult feces

From mid-May to early June, adults of tits and flycatchers breeding in artificial nest boxes were captured using the baffle method after the nestlings had hatched (Fan et al., 2021). The adults were placed in a cage (18 cm × 12 cm × 8 cm) lined with sterile craft paper by an experimenter wearing sterile latex gloves. These individuals usually defecated naturally within 5 min, and the feces were collected using sterile polyester swabs and placed in lyophilizing tubes (DNase/RNase-Free, Sterile) containing RNAlater (Servicebio, Wuhan, China). The lyophilizing tubes were gently inverted to completely immerse the feces in the RNAlater. The birds were released in place after they were ring banded. We obtained 49 fecal samples from adult tits across 32 nest boxes and 24 fecal samples from adult flycatchers across 19 nest boxes. The tubes containing feces were temporarily placed in a thermal bag with ice packs and then stored at −80°C until DNA extraction.

DNA extraction, PCR, and sequencing

Total DNA from collected fecal samples was extracted using a PowerSoil® DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA) according to the manufacturer's instructions. After DNA extraction, DNA integrity was tested using 1% agarose gels. We used the primers LCO-1490 (GGTCAACAAATCATAAAGATATTGG) and ZBJ-ArtR2c (WACTAATCAATTWCCAAATCCTCC) to amplify a 225-bp fragment of the cytochrome c oxidase subunit I (COI) barcode region (Folmer et al., 1994; Zeale et al., 2011). Sample-specific 8-bp barcodes were attached to the 5' ends of the primers. The PCR amplification was conducted in a final volume of 20 μl containing 4 μl of 5 × FastPfu Buffer, 2 μl of dNTPs (2.5 mM), 0.8 μl of primer F (5 μM), 0.8 μl of primer R (5 μM), 0.2 μL of Bovine Serum Albumin (BSA), 0.4 μl of FastPfu Polymerase, and 10 ng of DNA. The conditions for PCR were as follows: an initial denaturation step at 95°C for 3 min; 95°C for 30 s, 52°C for 30 s, and 72°C for 30 s, for a total of 45 cycles, followed by a final extension at 72°C for 10 min. Each PCR batch included 1–2 PCR blanks. Three separate PCR replicates were conducted for each extract, and these were then mixed after PCR amplification. After amplification, the PCR products were visualized by electrophoresis using 1% agarose gels and purified with Agencourt AMPure XP (Beckman Coulter, Indianapolis, IN, USA). The purified products were paired-end sequenced on an Illumina Miseq PE300 platform (Illumina Inc., SanDiego, CA, USA) by Beijing Ovison Gene Technology Co., Ltd (Beijing, China).

Sequencing data analysis

Raw sequences were split using QIIME v1.8.0 (Caporaso et al., 2010) based on the barcode tags. Paired-end reads were merged and quality-filtered by Pear v0.9.6 (Zhang et al., 2014). The sequences were processed using QIIME2 v 2020.6 (Bolyen et al., 2019). We removed duplicated sequences and chimeras using the Vsearch (Rognes et al., 2016) plugin in QIIME2, and the molecular operational taxonomic units (MOTUs) at a 97% similarity threshold were clustered (Vamos et al., 2017). We removed sequences that occurred only in a single sample and sequences with frequencies < 3. The MOTUs representing < 0.1% of the normalized sequences for each sample were removed to prevent the generation of potentially erroneous results (Bokulich et al., 2013). Representative sequences of each MOTU were compared with the reference sequences in the BOLD database (https://www.boldsystems.org/) and the Genbank database (https://www.ncbi.nlm.nih.gov/genbank/) to obtain taxonomic information. The taxonomic identification followed the criteria of Aizpurua et al., (2018) and Alberdi et al., (2018). Order-level taxonomic identifications were assigned at > 95% similarity values; family-level identifications were assigned at > 96.5%, and species-level identifications were assigned at 98% similarity. All of the identified species were checked manually. When one MOTU matched multiple species that shared the highest matching score, we downgraded the taxonomic resolution to the most common level. Those MOTUs not fulfilling the criteria or not matching any reference sequence were classified as unidentified.

Statistical analysis

For breeding time data, the first egg-laying date, incubation date, and hatching date were converted to Julian Days before analysis (1 April was defined as 1 each year (Fan et al., 2021)). For example, if the first egg-laying date was 1 May, this was recorded as 31. Generalized linear mixed models (GLMMs, Poisson distribution, log link, R package lme4) were used to examine differences in breeding time variables between tits and flycatchers, with species as a fixed variable. The first egg-laying date, incubation date, and hatching date were treated as response variables, and the year was treated as a random variable. We also calculated the breeding time overlap for each species. For example, the breeding time overlap for tits was calculated as the ratio of the co-breeding period with flycatchers to the total breeding period for tits. Due to the abnormal distribution of hatching rates (P < 0.05), we used the Mann-Whitney U test to assess the differences between tits and flycatchers.

For nest site data, Shapiro-Wilk tests were used to examine the normality of the distributions of nest site variables. If the characteristics were normally distributed, independent-sample t-tests were used to test for differences between tits and flycatchers. Otherwise, Mann-Whitney U tests were used (Appendix Table S1).

For diet data, we computed and plotted the rarefaction and extrapolation curves of prey species richness based on Hill numbers (q = 0) using the R package iNEXT (Hsieh et al., 2016), developed for presence/absence data. The 95% confidence intervals (CIs) were obtained by a bootstrap method based on 1,000 replications. The relative read abundance (RRA) and percent frequency of occurrence (%FOO) were used to quantify the prey composition at the order and family levels, respectively (Deagle et al., 2019). Kruskal-Wallis tests were used to test for differences in the RRA between the two species. Chi-square tests were used to test for differences in the %FOO between tits and flycatchers; arthropods with frequency > 2 were selected for the calculation. Independent-sample t-tests were used to analyze the differences in the individual diet species richness between the two bird species that were represented as the average number of arthropod species (MOTUs) consumed by each individual of each species. Nonmetric multidimensional scaling (NMDS) with Bray-Curtis distance was used to examine the degree of diet similarity between the two bird species, and an analysis of similarities (ANOSIM) test with Bray-Curtis distance was used to examine the difference in diet composition between the bird species using the R package vegan (Ramette, 2007). Finally, the diet niche breadth of each species was calculated using a standardized Levins index (BA) (Hurlbert, 1978).

To compare the niche occupied by the tit and flycatcher across multiple dimensions, we used a multivariate kernel density estimation method to measure interspecific niche overlap (Geange et al., 2011). Initially, we calculated niche overlap between the two species for each variable. Subsequently, we obtained the composite overlap across the three major ecological dimensions (breeding time, nest site, and diet) by averaging their corresponding variables. Finally, we obtained the overall niche overlap by averaging across all variables (Wang et al., 2021). We then used null model permutation tests to determine whether an observed overlap was significantly smaller than expected by chance. For each comparison, species labels were randomized 1000 times. The overlap between simulated kernel density distributions was compared with the observed overlap using t-tests. Significant differences indicated niche differentiation between the two species (Gotelli, 2000; Geange et al. 2011).

Except for the results of multivariate kernel density estimation (mean ± standard deviation), all other data were presented as mean ± standard error. All statistical analyses were conducted in R 4.1.2 (R Core Team 2022). The significance level was set at 0.05.