Reassessment of French breeding bird population sizes using citizen science and accounting for species detectability

Sampling protocol

EPOC-ODF (French structured estimation of breeding bird population size) is a French CS monitoring scheme based upon 5-min point counts, where observers are tasked to point locations of recorded individuals, either through visual or auditory detection. Birders can register their field observations directly using the NaturaList smartphone application or transcript later on the data portal Faune-France (www.faune-france.org). The survey locations corresponded to the centroids of a 2 × 2 km grid, selected from a random sampling. Each location has to be visited three times during the breeding season, from March to June, each consisting of three successive 5-min point counts, to limit chances of duplicated counts while being less demanding in observation effort (Fuller & Langslow, 1984). After completion, i.e., nine visits during a breeding season, surveyed sites are removed from the sampling pool for the subsequent year, to maximise the number of sites surveyed. See Materials S1 for more details about the sampling design.

Over the 2021 and 2023 breeding seasons, 276 distinct species were encountered over 27,156 complete checklists collected over 3,873 pre-selected locations (Fig. 1) by 520 observers. Sampling effort is monitored through local associations tasked to recruit volunteers. The primary focus of the scheme being the monitoring of common breeding bird species, we decided to constrain the number of species considered viable targets of this scheme to 103 out of the 276 species contacted. We narrowed our study to 63 species out of the initial set of 103, comprising only those recorded at a minimum of 150 distinct locations (3.9% of total locations), to have a sample size allowing to reach model convergence. We also applied a temporal filter that considered both observed activity during the breeding season and expert opinion to capture the breeding phenology of each targeted species and exclude possible early or late migrants from population size estimates (see Materials S2.1).

PCA reduction of environmental covariates

For bioclimatic data, we used 19 variables from WorldClim at 1 km resolution (Fick & Hijmans, 2017), on which we applied a Principal Component Analysis (PCA), keeping the first three axes (82.3% of explained variance), to limit multicollinearity through orthogonal transformation (Cruz-Cárdenas et al., 2014).

We used habitat cover data from Theia OSO at 10 m resolution (Thierion, Vincent & Valero, 2022) and aggregated it according to two different scales: (1) a seven-class corresponding to habitat type (Urban, Annual crops, Perennial crops, Pastures, Grasslands, Forests, Water body/Mineral surfaces) and (2) three-class (Open, Forests, and Artificial) in regards to overall effect on detectability (Fig. 2). Additionally, we conducted PCA dimension reduction on the seven-class aggregation, retaining three of the six PCA (54.71% of explained variance) axes depicting environmental gradients for (i) forest-to-open-field cultures; (ii) open-field cultures-to-pastures and (iii) perennial crop-to-urban Materials S3 for the workflow pipeline and habitat cover aggregation. Distances to roads were measured from ROUTE 500 (Cote et al., 2021). Environmental covariates were extracted over a 500 m buffer radii upon registered observer location (Fig. 2). These distances were chosen according to mean dispersal distances and home range sizes in common European birds (Paradis et al., 1998). The three-class habitat covers were collected upon 100 m circles radii to assess immediate habitat types that could hinder species detection. Whenever the exact location was unavailable, we used the centroid of sightings as a proxy for observer location (Materials S4). We used environmental data collected from a prediction grid covering France at a resolution of 2 × 2 km for PCA dimension reduction. Outcomes from this initial PCA were used to transform environmental data collected from surveyed locations through PCA projections.

Modelling framework

We used Hierarchical Distance Sampling (HDS) models to estimate the abundance of the target species while accounting for uncertainty arising from the observation process (Chandler, Royle & King, 2011; Kéry & Royle, 2015). We applied a right-side truncation of 5%, removing observation distances above the 95% quantile, for each targeted species to remove extreme distance values for model robustness (Buckland et al., 2001). Then, we divided observation distances into five proportional bin classes based on the maximal observed distance. Models calibration and assessment were done using unmarked 1.2.5 R package (Fiske & Chandler, 2011). Effort covariates were accounted for by incorporating the Julian date and the hours of list realisation (as minutes from sunrise), see Table 1.

Distance sampling key functions, depicting detection probabilities fall-off given distance of observation (Buckland et al., 2001), were chosen between half-normal and hazard-rate based on Akaike Information Criterion (AIC, Akaike, 1974), with other states kept constant.

We based our modelling framework on a secondary candidate set strategy (Fig. 3), where the detection and availability states of our HDS were fit according to the set of the first candidates while others were kept constant (Morin et al., 2020). For the Poisson process underlying abundance distribution, we used a single model consisting of retained covariate PCAs axes (Table 1). See Table S5.1 for the number of times where each sub-process was included in the final candidate sets.

HDS population size estimates were obtained by averaging retained secondary candidate sets models, based on their relative model performance using AICc (Fig. 3). We excluded the Eurasian Sparrowhawk (Accipiter nisus), the Meadow Pipit (Anthus pratensis) and the Coal tit (Periparus ater), from model averaging and exclusively relied upon prediction from best final models owing to substantial differences observed among their secondary candidate sets models.

Model goodness-of-fit was assessed using an oveardispersion coefficient metric ( $\hat{C}$; Johnson, Laake & Ver Hoef, 2010). We used the chi-square metric as the discrepancy measure between observed and expected counts. Computed $\hat{C}$ corresponds to the ratio between the chi-square obtained from the fitted model to the mean of bootstrapped chi-squares obtained from simulated datasets based upon estimated parameters (Kéry & Royle, 2015). All models were fit according to a Poisson (P) distribution after top model assessment and calculation of $\hat{C}$, secondary candidate sets with $\hat{C}$ top models exceeding 1.2 were calibrated using a negative binomial (NB) distribution (Payne et al., 2018). For a global overview of our modelling approach, see Materials S5. Out of 63 species, we excluded nine species from the analysis; three exhibited signs of underdispersion with $\hat{C}$ values less than 0.9 while six had $\hat{C}$ values exceeding 1.5 (Payne et al., 2018), showing signs of overdispersion, despite being calibrated using a negative binomial distribution, see Table S2.1 for more details.

We assessed the robustness of our estimations to the exclusion of one year of data, corresponding to a third of the global dataset. We compared population size estimates from EPOC-ODF data collected over 2021–2023 to estimates obtained from EPOC-ODF data collected over 2021–2022. Using the 2021–2022 subset, we estimated the population sizes of 30 species, detected in at least 150 distinct sites. From these 30 species, seven mean population sizes estimated using 2021 to 2023 data were outside the confidence intervals estimated from 2021 to 2022 data, with a slightly smaller population size estimated (Table S2.2) overall highlighting robust estimations.

Trimming of HDS population size estimate: assessment of model extrapolation

Population sizes were obtained by summing predicted values over the prediction grid. As we intend to predict over a large surface, novel environmental conditions may arise, leading to possible dissimilarities between the environmental gradient collected at survey sites and the environmental gradient over novel conditions (Elith & Leathwick, 2009). Mesgaran, Cousens & Webber (2014) described two types of extrapolation; (i) novelty type I (NT1) where projected points (i.e., prediction grid) are outside the range of individual covariates collected by the sampling scheme and (ii) novelty type II (NT2) depicting the case when projected points are within univariate range but constitute novel combinations between covariates.

We trimmed predicted values (Fig. 4) over prediction grid cells showing signs of NT1 extrapolation, to a threshold value determined using a Tukey fence (Tukey, 1977), estimated from the distribution of predicted values with k = 1.5. Extrapolation assessments were done using dsmextra 1.1.5 R package (Bouchet et al., 2020). We used this post-prediction treatment to assess population size estimates stability. We measured the coefficient of variation, corresponding to the ratio between the standard deviation and the mean of the “untrimmed” and “outlier-trimmed” estimated range uncertainty. Large coefficients of variation imply great discrepancies in confidence intervals of untrimmed and outlier-trimmed estimated uncertainty intervals. This is mainly caused by the spatial filtering from the extrapolation assessment highlighting smaller geographic regions with similar environmental conditions from sampled ones and the trimming of predicted abundance outliers. Species with a coefficient of variation exceeding 30% were removed from the comparison of ArGeom and HDS population size estimates (Table S2.3).

Comparison of ArGeom and HDS estimated population sizes

For comparable estimates between ArGeom and HDS approaches, we restricted the prediction grid area species-wise for HDS estimation according to the distribution of their known breeding locations, collected over a 10 × 10 km grid during the previous French atlas (Issa & Muller, 2015). To estimate breeding populations of species for which male identification was possible, either male vocalisations or visual distinctions because of sexual dimorphism, an ad-hoc filter was applied (Table S2.1), resulting in HDS estimates reflecting the male counts for those species.

As the ArGeom approach estimated species bird population sizes as a number of breeding pairs (Roché, Muller & Siblet, 2013), for species where male identification in the field was impossible (no sexual dimorphism), we used all available data, after applying the phenological filter, and divided HDS estimates by two for comparable estimates with ArGeom population sizes. After retrieval of ArGeom estimates from the previous atlas (Issa & Muller, 2015), we updated these estimates using recent population trend estimates derived from the French breeding bird survey (FBBS; Jiguet et al., 2012) data spanning 2012–2023 (Table S2.1). Given the absence of a mean estimate in the ArGeom approach, we approximated it using the midpoint between the maximum and minimum estimated (Fig. 5).

To study the differences between the two approaches, we measured ${\delta }_{mean}$ corresponding to the percentage of the difference between HDS and ArGeom estimates.

$${\delta }_{mean}={\displaystyle \frac{\left(Estimate{s}_{ArGeom}-Estimate{s}_{HDS}\right)}{\left(Estimate{s}_{ArGeom}+Estimate{s}_{HDS}\right)/2}.}$$

Study of variation of estimated population sizes between the two approaches

As species detectability stems from physical traits and vocalisations, phylogenetic related species tend to have the same detectability (Johnston et al., 2014; Sólymos et al., 2018). We calibrated a Phylogenetic Generalised Linear Mixed Model (PGLMM) using the phyloglmm 1.0 (Li & Bolker, 2019) R package. We study ${\delta }_{mean}$ variations across species while implementing a random effect covariance structured based on phylogenetic relatedness using phylogenetic distances retrieved from Burleigh, Kimball & Braun (2015). The PGLMM model was calibrated using (i) extracted detection probabilities from the availability state estimated through HDS (Fig. 3 and Materials S5) after model averaging of the final candidate sets models in regards to AICc scores, and (ii) ArGeom uncertainty as fixed variables. For ArGeom uncertainty, corresponding to the difference between maximal and minimal estimated values, we relied on the decimal logarithm to limit variation in ${\delta }_{mean}$ solely due to different population size magnitudes.

Response weights consisted of normalised weights from the inverse of uncertainty around FBBS trends between 2012 and 2023 (Table S2.1), divided by the mean to limit excessive weight attribution and facilitate model convergence.

Species trends over 2012–2023

From 2012 to 2023, out of 63 bird species, 15 showed a significant decrease ($\overline{x}$ = −22.79% $\pm$ 14.84) in total population size, while 16 showed a significant increase ($\overline{x}$ = 28.02% $\pm$ 22.52; see Materials S2 for species-related FBBS trends).

HDS population size estimations

Out of the 54 species with acceptable values of overdispersion ( $\hat{C}$) using the HDS approach, we excluded eight species showing large discrepancies in population size estimates (Materials S2, Table S2.3) between pre- and post-prediction treatment (Fig. 2.).

Out of the remaining 46 species used for comparison between ArGeom and HDS estimates, HDS models showed acceptable values of overdispersion ( $\hat{C}$) ranging from 0.94 to 1.2 ( $\overline{x}$ = 1.07 $\pm$ 0.06) for 38 species calibrated using a Poisson distribution and 1.09 to 1.47 ( $\overline{x}$ = 1.27 $\pm$ 0.13) for eight species calibrated using a Negative binomial distribution.

Population size comparison between ArGeom and HDS

Across all species, estimated mean density ranges from 0.09 to 27.51 individuals per square kilometre, while ArGeom range uncertainty varies from 3.9 to 6.69 on the decimal logarithm scale corresponding to variations from 7,920 to 4,850,000 in estimated number of pairs. See Tables S6.1–S6.3 for more details about species estimated population size according to ArGeom and HDS approaches.

A comparison between updated ArGeom and HDS estimated population sizes showed that HDS estimates were higher than ArGeom for 30 of the 46 species tested (Table S6.1). Our results suggest lower estimates from ArGeom ( ${\delta }_{mean}$ < −0.2), either for open habitat specialists such as European Stonechat (Saxicola rubicola), European Goldfinch (Carduelis carduelis) or Eurasian Linnet (Carduelis cannabina) than for forest generalists such as Great Spotted Woodpecker (Dendrocopos major), Blue Tit (Cyanistes caeruleus) or Eurasian Blackcap (Fig. 6A) with ${\delta }_{mean}$= −0.629 $\pm$ 0.4, over 22 species. Species whose estimations were similar ( ${\delta }_{mean}$ ∈ [−0.2;0.2], with on average ${\delta }_{mean}$ = −0.006 $\pm$ 0.096, over 17 species) between ArGeom and HDS included species such as Eurasian Wren (Troglodytes troglodytes), Great Tit (Parus major), European Robin (Erithacus rubecula) and European Nuthatch (Sitta europea). Fewer species, mainly characterised by greater maximal observation distances, such as Common Cuckoo (Cuculus canorus), Eurasian Hoopoe (Fig. 6B) and European Blackbird (Turdus merula) had higher population sizes estimated by ArGeom approach compared to HDS ( ${\delta }_{mean}$ > 0.2, with on average ${\delta }_{mean}$ = 0.41 $\pm$ 0.148, over seven species; see Materials S6 for population size comparison table).

Results from the PGLMM (Fig. 7A) showed an overall significantly lower ArGeom population size estimates (−0.209, with 95% CI [−0.321 to −0.097], Pval = 0.001), as well as a significant positive effect of ArGeom range intervals (0.167, with 95% CI [0.059–0.276], Pval = 0.003) on the differences between the two approaches. Species detection probabilities had no significant effect (0.097, with 95% CI [−0.043 to 0.237], Pval = 0.176) on ${\delta }_{mean}$ variation.

Marginal effect plots from the PGLMM model showed that the mean response of ${\delta }_{mean}$ over species detection probability was predominantly negative, ranging from −0.45 to −0.1 (Fig. 7B), for ArGeom uncertainty. This showed that ${\delta }_{mean}$ tended towards the convergence of population size estimates (Fig. 7C) for species with larger estimated interval ranges. There were no signs of multicollinearity (VIF < 5; James et al., 2013) between the two variables.

Potential consequences for community-level assessments

A recent study about long-term effects of climate and land use changes on bird communities (Gaüzère et al., 2020) showed that both generalist cold-dwelling species, such as the Common Chiffchaff (Phylloscopus collybita) or the Eurasian Blackcap, and warm-dwelling species, such as the Common Nightingale (Luscinia megarhynchos) had the most substantial negative and positive contributions to the trend in Community Thermal Index (CTI), a community-weighted index representing the realised thermal niche of a community based upon species relative abundance and species thermal indices (STI). In the present study, these species tended to have lower population sizes estimated when the detection process was omitted compared to estimates based on our modelling approach. As a result, this could affect the estimations of their contribution to the calculation of community-weighted mean indices, such as CTI, and therefore bias the estimation of the trend in community thermal response and subsequent studies of aggregated indices, which are known to display large regional variation (Rigal & Knape, 2024). We, therefore, suggest that considering the detection process in studies relying on community-weighted indexes by species’ relative abundances could be as important as it is for estimating population sizes.

Community indices such as species diversity (Ricotta, 2005) and functional diversity (Villéger, Mason & Mouillot, 2008; Gaüzère et al., 2019) are commonly use species relative abundance as a basis, without taking into account the detection process (Pillar & Duarte, 2010), despite multiple studies showing it could affect community indices inference (Tingley & Beissinger, 2013; McNew & Handel, 2015; Jarzyna & Jetz, 2016; Richter et al., 2021).

Conservation implications

Our study also suggested that lower or higher population sizes estimated from ArGeom were not randomly distributed among species according to their conservation status. Out of the 46 species estimations used in the comparison analysis, 10 had an unfavourable conservation status in France (i.e., lower than Least Concern, LC; UICN France et al., 2016).

Among these species of conservation concern, two species, European Greenfinch (Carduelis chloris) and European Turtle Dove (Streptoptelia turtur), showed no signs of difference in their population sizes. By contrast, five species, European Stonechat, Barn Swallow (Hirundo rustica), Red-backed Shrike (Lanius collurio), Eurasian Kestrel (Falco tinnunculus) and Willow Warbler (Phylloscopus trochillus) considered as NT (Near Threatened) and three species, Eurasian Linnet, European Goldfinch and European Serin (Serinus serinus) considered as VU (Vulnerable) had lower population sizes estimated from ArGeom than from HDS approach (NT: ${\overline{\delta }}_{mean}$ = −0.608 $\pm$ 0.217 and VU: ${\overline{\delta }}_{mean}$ = −0.667 $\pm$ 0.146). Our results showed that these species may need a reevaluation of their conservation status and highlight the need to rely on hierarchical models taking account of the detection process in ecological inferences, given that potential misclassification of population conservation status may arise from process noise and observation error (Connors et al., 2014). As conservation policy decisions depend on uncertainty levels (Williams, 2003; Freckleton, 2020), assessing measurement error through the integration of the detection process (Nichols et al., 2011) could provide more reliable ecological inferences (Guillera-Arroita et al., 2014). CS schemes are becoming more and more a reliable source of data to ensure biodiversity monitoring (Chandler et al., 2017) and can, through standardisation (Buckland & Johnston, 2017; Johnston et al., 2019), contribute to the calibration of data-hungry models such as hierarchical models for reliable ecological inferences (Isaac et al., 2020; Kéry & Royle, 2021; Johnston, Matechou & Dennis, 2022).

Comparison to other European countries

Another way to assess the relevance of the two estimation approaches would be to compare their population size estimates to the ones obtained from other European countries, using a ratio between countries to produce comparable estimates. Such an approach should however be used with caution because it would be limited by comparability in habitat repartitions or biogeographical considerations among different European countries. To go further into inter-country comparisons, we relied on the German population sizes estimated for the previous European Bird Directive (BirdLife International, 2021) obtained from both point count and territory mapping methods (Gedeon et al., 2015). For abundant species such as the Blackcap ( ${\delta }_{mean}$ = −0.29; German population size expressed in millions of pairs = [7.17–9.49]), both approaches led to similar results than German population estimates, while HDS estimates were closer to German population sizes for the Firecrest (Regulus ignicapilla; ${\delta }_{mean}$ = −0.46; [1.92–2.85]) and the Blue Tit ( ${\delta }_{mean}$ = −0.42; [5.01–7.41]). For species with higher population sizes estimated by ArGeom than HDS ( ${\delta }_{mean}$ > 0.2), the Common Cuckoo ( ${\delta }_{mean}$ = 0.25; [0.58–0.95]) and the Corn Bunting (Emberiza calandra; ${\delta }_{mean}$ = 0.6; [0.25–0.44]) showed estimates of German populations closer to the HDS than the ArGeom approach. Finally, for the Common Whitethroat (Curruca communis; ${\delta }_{mean}$ = 0.42; [0.93–1.47]), the German population size is closer to ArGeom estimates (see Materials S6.4 for additional information).

Regarding magnitudes, both approaches produced similar estimates compared to German ones. However, due to different sampling and modelling methods used, these formal comparisons, although informative, need to be more fully satisfying and highlight the discrepancies in sampling and analytical methods across the European continent (Keller et al., 2020). Such differences could be accounted for, either by (i) a global standardisation of schemes, as promoted by the PECBMS (Pan-European Common Bird Monitoring Scheme; Brlík et al., 2021) for species trends, but also (ii) through the use of Integrated Models (IM) capable of mobilising data from multiple and somewhat heterogeneous sources (Isaac et al., 2020; Zipkin et al., 2021).

Study limitations

Our approach relies on data collected from the EPOC-ODF structured CS schemes, providing data with repeated visits. However, as is, the frequentist framework of unmarked R package (Kellner et al., 2023) does not permit inferences on social species occurring in large flocks. Taking account of social species during the breeding season (corresponding to 1/10th of the scheme targeted species) would therefore require a Bayesian framework to include the effect of flock size on species detectability (Clement, Converse & Royle, 2017).

Given the timeframe and the sampling design, i.e., all sites are not visited every year to maximise the number of total surveyed locations, it is not possible to estimate species demographic parameters, such as survival and recruitment (Sollmann et al., 2015; Schmidt & Rattenbury, 2018). We also assumed a sex ratio of 1:1 for species without sexual dimorphism, during the breeding season, which could potentially bias estimates for species deviating from this assumption. Taking account of species population structure requires frameworks such as Integrated Population Models (IPM; Schaub & Ullrich, 2021) and specific data collection (King, 2014), for instance, bagging or nest surveillance.

As obtaining relevant predictions of species abundance over unsampled environmental conditions was one of our main methodological challenges, we used environmental data condensing habitat information (Tredennick et al., 2021). To fit our statistical framework, we assumed that most bird species would interact with their habitat following a linear relationship (see Fig. 3). We therefore used PCA reduction to summarise species linear responses to national-scale habitat gradients including forest-to-open-field cultures, open-field cultures-to-pastures and perennial crop-to-urban habitats (see Materials S3). PCA reduction permits model convergence by condensing complex habitat structures to a small number of environmental covariates, though it could bias estimates of species thriving in a specific habitat restricted to the extreme edge of the sampled gradients. Other methods such as Spatially Varying Covariates models (SVC; Gelfand et al., 2003) could be used to better account for habitat structure complexity across spatial gradients (Thorson et al., 2023).

Previous studies have shown that unaccounted variations in species availability, considering a constant detection probability or unmodelled variations, could lead to substantial bias in estimated abundance (Link et al., 2018; Barker et al., 2018; Duarte, Adams & Peterson, 2018). N-mixture biassed estimations can be linked to non-assessment of the sampled area, where a smaller or greater sampled effective area could lead to under- or overestimation (Kéry & Royle, 2015). In our study, as we relied on distance sampling methods, we define an effective sampled area, based upon collected observation distance, but we also assumed that individuals considered exposed to the sampling (i.e., ‘statistically’ available for modelling) could still be undetected due to small species home ranges or plot-specific habitat cover (Chandler, Royle & King, 2011; see Table 1 and Fig. 2, for covariates used to model species detectability and Materials S5 for model formulation). Despite such consideration, for the HDS model, we assumed that detected individuals were homogeneously distributed over the sampled area. Violating this assumption could lead to within-sample variation that needs to be accounted for, otherwise leading to biassed estimates (Mizel, Schmidt & Lindberg, 2018).

Another potential drawback relies on the quantity of data collected through this structured CS scheme. Over the same breeding season, the semistructured scheme EPOC without temporal replicates nor fixed location requirements collected three times the amount of complete checklists as the structured EPOC-ODF scheme, highlighting CS trade-off of scheme standardisation upon data collection over spatial and temporal scales (Devictor, Whittaker & Beltrame, 2010). One way to address this trade-off would be to apply data integration methods mobilising multiple data sources to be used for ecological inferences (Zipkin, Inouye & Beissinger, 2019; Zipkin et al., 2021), either by estimating abundance of less recorded species through trait-based associations (Callaghan, Nakagawa & Cornwell, 2021, 2022; Robinson et al., 2022) or by constructing joint likelihood functions (Fithian et al., 2015; Fletcher et al., 2019).