Data Collection

We built an IPM and fitted parameters with 4 datasets: (1) population survey data, (2) nest survey data, (3) duckling survival data, and (4) adult band-recovery data. Our study area covers the entire country of the Netherlands, as all data were collected as part of national monitoring programs or citizen science projects. We limited our analyses to the years 2003–2020 due to a limited sampling effort in the years before 2003 for most datasets.

Population survey data. Breeding bird populations have been monitored in the Netherlands as part of a national Breeding bird Monitoring Program (BMP) starting in 1984 (Boele et al. 2014). BMP is based on intensive territory mapping in fixed study plots carried out by well-trained volunteers and professionals who follow a standardized protocol. Territory mapping is based on a large, and annually constant, number of field visits (5–10 between March and July depending on species) in which all birds with territorial behavior (e.g., pair bond, display, nests) are recorded on maps. Species-specific interpretation criteria are used to determine the number and exact locations of “territories” at the end of the season. The counts of 1990 constitute the index baseline. Index values for the following years are then calculated as the relative difference between that year's counts and those of the reference year (Sovon 2019).

Nest survey data. Nest data for Dutch Mallards were retrieved from the National Nest Record Scheme (Bijlsma et al. 2020). Mallard nests are mainly reported as “bycatch” during nest surveys of grassland-breeding waders. To calculate clutch size and egg hatch rate, incomplete clutches were excluded by using only successfully hatched nests. In total, there were 4,415 observations of clutch size with a range of 149–402 nests observed annually. For egg hatch rate per successful clutch, we calculated for each year the total number of eggs laid and the number of eggs hatched in nests where at least one egg had hatched. In total 2,383 nest observations were available, with the number of observations ranging between 54 and 282 nests per year.

Duckling survival data. We used data for Mallard females with broods collected through a citizen science project to estimate annual duckling survival rates. Contributors to this project were asked to report brood size and brood age on a mobile application specifically developed for this purpose. From 2018 onward, contributors were asked to make repeated observations of the same broods to allow for the calculation of duckling survival. We therefore used data for the years 2018–2020, as only for these years were there enough repeated observations of broods to estimate duckling survival. The final dataset included 2,825 observations of 1,212 broods distributed over the three years.

Adult band-recovery data. Annual adult survival was estimated using data for Mallards that were banded between 2003 and 2020 and were later recovered and reported dead. In the Netherlands, Mallards are banded mostly in traditional duck decoys (Karelse 1994). Since the year 2000, 20,133 Mallards have been banded, of which 1,233 were recovered dead (Buijs and Thomson, 2001, van Noordwijk et al. 2003). We only included Mallards that were banded in the breeding season (March–July) to exclude winter migrants that breed in other countries. This dataset included 2,615 Mallards, of which 203 were recovered dead.

Model Structure

The IPM was based on a prebreeding census life cycle model that distinguishes 2 life-stages: 1-year-old females and females aged 2 years and older. Male Mallards were excluded from the model as they outnumber the females in this population (Eygenraam 1957, Boele et al. 2022) and, therefore, do not pose limits on the production of offspring. The cause of this sex ratio imbalance remains unclear (Boele et al. 2022), and because previous studies of European Mallards have reported only minor sex- and age-related differences on survival (Gunnarsson et al. 2008, 2012), we did not use sex- or age-specific vital rates. We therefore assumed that all individuals contribute equally to the population of breeding females in the following year.

Likelihoods for the Datasets

We first describe how the population survey data are linked to annual reproductive and survival rates, whereafter we explain how these are modeled using the individual survey datasets.

Population survey data. We used a stage-structured state-space model (SSM) to formulate the likelihood of the population census data. The SSM links the demographic rates of reproduction and survival to the underlying true population size (N_t) and the observed population index of year t, where t = 1 corresponds to the values for the year 2003. In this model, the population size of a given year N_t+1 is determined by 2 inputs: (1) the number of female Mallards that were produced in the preceding breeding season and have survived until the current one (N_R,t+1), and (2) the number of adult females that survived the interval between these two periods (N_S,t+1). This state process can be described with:

Here, ρ_t is the annual reproductive rate, defined as the number of females born per adult female that survive until the end of a breeding season, while the annual probability that an individual survives from one breeding season to the next is given by φ_t. Note that Eqs. (1) and (2) use the normal approximations of respectively the Poisson and binomial distribution, because the population survey data consist of continuous index values instead of discrete counts (e.g., Besbeas et al. 2002, Brooks et al. 2004). Annual reproductive rate is calculated as the product of 4 lower vital rates, each of which is conditional on the one before it. These are (1) nest success (π_t): the probability that a female successfully hatches at least one egg during the breeding season; (2) clutch size (C_t/2): the number of eggs laid per nest, divided by two to represent the female eggs only; (3) egg hatch rate (E): the probability that an egg hatches; and (4) duckling survival (D_t): the probability that a hatched duckling survives until fledging age. Next, the observation process of the SSM relates the true underlying index values, as determined by the state process, to the values observed by the data with:

We here denote the annual observed population index with values y_t where y₁,..., y₂₀ and assume that each observation y_t is made with constant variance σ_y². The population growth rate between breeding seasons was calculated for the intervals 2003–2004 to 2018–2019 as:

Nest survey data. The annual probability that a Mallard successfully completed a single nest (π_1,t) was modeled using the Mayfield estimate of daily nest survival that relates the number of nest failures (y_fail,t) to the number of days those nests were exposed to potential failure (y_expo,t; Mayfield 1975). The daily failure rate of nests (DFR_t) was assumed to be constant during the incubation period and modeled using the following binomial:

Parameter π_1,t was then calculated by raising (1 – DFR_t) to the appropriate power that represents the time interval in days between nest initiation (first egg laid) and nest completion (last egg laid; Mayfield 1975). For Mallards this corresponds to a nesting period of ∼37 days. Because Mallards can attempt multiple nests after initial nest failure (Arnold et al. 2010), true annual nest success is underestimated by π_1,t, and to account for this, we used the formula provided by Hoekman et al. (2002) to derive true annual nest success (π_t):

Here, p₁ is equal to breeding incidence while p₂ through p₆ are described by p_z = max{1 – 0.26(z – 1), 0}, where z is the z^th nest attempt (Hoekman et al. 2002). This formula describes the cumulative nest success probability of a female during the breeding season. Observations of clutch size in individual Mallard nests (y_c,i,t) were modeled using a Poisson distribution with mean annual clutch size (C_t):

Egg hatch rate E was assumed to be invariable throughout the years and modeled per nest i as the number of eggs hatched (y_eh) given a number of eggs laid (y_el) using a binomial distribution:

Duckling survival data. For the years 2018–2020 we modeled duckling survival from hatching to eight weeks of age using repeated observations of Mallard brood counts. Ducklings fledge approximately 56 days after hatching and around this time, daily survival rates approach 1.0. Therefore, we assumed that mortality of juveniles post-fledging until the end of the breeding season was negligible. Observations of brood size (y_b) from a particular female at day τ₂ were modeled as a Poisson distribution with mean brood size reported during the previous observation of the female at day τ₁. This was then multiplied by the daily duckling survival rate (φ_d) and raised to the power of the time interval between the two observations (τ₂ – τ₁):

Duckling survival is known to strongly depend on duckling age and to be highly variable between years. To confirm this, we performed model selection on the duckling survival model and used the Deviance Information Criterion to compare the age- and year-dependent model to a null model and models only incorporating one of the two variables ( Supplementary Material Table S4 (duac020_suppl_supplementary_material.docx)). As the model with both variables had the most support, we modeled duckling survival as:

Variable Age is here the average of the estimated duckling ages at days τ₁ and τ₂, and to account for annual differences in survival, each year was given its own intercept (β_D0,_t) and slope (β_D,t). The probability that a duckling survives the first eight weeks of its life (D_t) was then calculated as the product of the survival probabilities for days 1 to 56, which are given by the inverse logit functions of Eq. (3) for Age = 1,..., 56:

Adult band-recovery data. Adult survival of Mallards was modeled according to the dead recovery model introduced by Catchpole et al. (1999) using code provided by Brooks et al. (2004). Like van den Bremer et al. (2015), we considered the breeding season the census period. The likelihood of the band recovery data is determined by two parameters. These are the probability that an individual survives a given year (φ_t), and the probability that an individual is recovered dead and reported as such in a given year t (r_t). These parameters combined to produce a multinomial distribution that denotes the probability of recovering an individual dead in year j that was previously released in year i. This probability model was then linked to the observed recovery data. These were presented in the upper triangular matrix m_i,j, which lists for each released group of Mallards i the expected number of dead recoveries in the following years j. The number of banded Mallards that were never recovered dead are given in the final column j = T. The likelihood of observing m_i,j was then optimized for parameters φ_t and r_t according to:

where , representing the expected number of banded birds that are never seen again during the time of the study. For more details on this model see Catchpole et al. (1999).

Priors and Model Specification

We adopted a hierarchical approach similar to Cleasby et al. (2017) to get estimates of annual variation in DFR_t, C_t, φ_t, and r_t. Here, annual estimates of demographic processes were assumed to be drawn from normal distributions for each parameter, with means and variances denoted by respectively β_{vital rate, 0} and σ²_{vital rate}. We additionally incorporated a linear trend for DFR_t and C_t described by the coefficients β_{vital rate, 1} multiplied by Time_t, which represents the standardized years. We did not add linear trends for φ and r, as external model selection in program MARK (White and Burnham 1999) did not support time-dependent adult survival. Similarly, we did not estimate separate survival probabilities for Mallards of different ages due to negligible differences between age classes ( Supplementary Material Table S3 (duac020_suppl_supplementary_material.docx)). This yielded the following formulas:

Here, , , , and . We selected uninformative priors for all model parameters with σ_r ∼U(0,10), σ_φ ∼U(0,10), σ_y ∼U(0,10), and σ_C ∼U(0,10). For both the intercepts and regression coefficients for the hierarchical modeling we used the priors β_{vital rate} ∼ N(0,100). Because there were no data for duckling survival for the years 2003–2017 we used the vague priors D_t ∼U(0,1) for t = 1,...,15. For 2018–2020 with available data we used βD,t∼ N(0, 100) and βD₀,t ∼ N(0, 100) to model duckling survival as a function of duckling age (Eq. 3). Priors for the egg hatch rate and the latent index value for the first year were selected as E ∼U(0,1) and N₁ ∼ N(101,1). We constructed the IPM using the program JAGS (Plummer 2003) and ran it with the r2jags package (Su and Yajima 2015) written for the programming language R (R Core Team 2020). We maximized the joint likelihood of the model using Monte Carlo Markov Chain simulation using 3 chains, 200,000 iterations, a burn-in of 100,000, and a thinning interval of 100. We checked for parameter estimate convergence (R^{^} < 1.02) using the Gelman-Rubic convergence statistic (Brooks and Gelman 1998). To visualize the effect of the vital rates on the projected population growth rate, we used the R package popbio (Stubben and Milligan 2007) to estimate λ for a range of values for ρ and φ.

Sensitivity and LTRE Analyses

We performed life-stage sensitivity analysis as described by Wisdom et al. (2000) to estimate the sensitivity of λ to the lower-level vital rates. First, we used the IPM estimates to generate probability distributions for the 5 vital rates and their annual variabilities. To this end we used β-distributions based on vital rate variances for all parameters except clutch size, for which we used a uniform distribution based on the range of observed annual values. We then sampled each distribution to create 10,000 replicates of the 5 vital rates. Finally, we used the vitalsens function in the popbio package to calculate the sensitivity of λ to the lower-level vital rates across the range of realistic parameter values. For more details on this approach see Wisdom et al. (2000).

To get an indication of which vital rates have changed in recent decades, and thus may have contributed to the population decline, we searched for historic values in the literature. We then used the Life Table Response Experiment (LTRE) approach, by multiplying the difference between current and historic values of vital rates with their λ-sensitivity value, to estimate how much their change could have contributed to a decline in population growth (λ) over time (Horvitz et al. 1997, de Vries et al. 2022). R code and data to reproduce our analyses can be found in Wiegers et al. (2022a, b).

Estimates of Demographic Rates

The model predicted a mean population growth rate of 0.98 (95% credible interval [CI]: 0.97–1.00), corresponding to an annual population decline of 2% (Figure 2 and  Supplementary Material Table S1 (duac020_suppl_supplementary_material.docx)). Predicted values for the unobserved true population index N_t conformed well with the observed annual index values ( Supplementary Material Figure S1 (duac020_suppl_supplementary_material.docx)). Estimates of annual clutch size were relatively invariable (Figure 3A) while estimates of mean nest success fluctuated substantially throughout the study period (Figure 3B). We found no significant temporal trends for either clutch size or nest success ( Supplementary Material Table S1 (duac020_suppl_supplementary_material.docx)). The estimates of duckling survival for 2018, 2019, and 2020, the only years for which empirical data were available, were 0.25 (95% CI: 0.18–0.32), 0.16 (95% CI: 0.12–0.21), and 0.27 (95% CI: 0.24–0.30), respectively. Duckling mortality decreased rapidly with duckling age and became almost negligible at four weeks of age ( Supplementary Material Figure S2 (duac020_suppl_supplementary_material.docx)). Mean duckling survival during the entire study period was within the observed annual variation in the years for which data were available. Model estimates for annual duckling survival and reproductive rates had large uncertainty margins for the years 2000–2017 due to uncertainty in duckling survival for those years (Figure 3C). Lastly, adult survival did not vary substantially among years (Figure 3D).

FIGURE 2.

A contour plot depicting how annual adult survival and reproductive rate determine population growth rate for the Dutch population of the Mallard. Black dots represent values as estimated by an IPM for the years 2003–2020. The white square indicates the mean population growth for the entire period (λ = 0.98). Contour lines represent annual projected population growth rate values.

Effect of Demographic Rates on Population Growth

Reproductive rate was strongly correlated with duckling survival (r = 0.64, df = 16, P = 0.004, Figure S3A), but not with clutch size (r = 0.009, df = 16, P = 0.97; Figure S3B) or nest success (r = 0.04, df = 16, P = 0.89; Figure S3C). The sensitivity analysis showed that the population growth rate was the most sensitive to duckling survival. Additionally, duckling survival had the largest LTRE contribution that far exceeded those of nest success and adult survival (Table 1), meaning that the drop in duckling survival (from an estimated 0.50 in the 1950s to 0.18 in the current analyses) was estimated to have the largest effect on λ.

Vital Rate Estimates

Duckling survival. Annual estimates for duckling survival as observed in 2018–2020 were relatively low compared to both historical estimates from the Netherlands and studies in other regions. It should be noted that, in reality, duckling survival was probably even lower than the estimated 16–27%, as observers were unable to monitor broods from the exact moment of hatching. Duckling mortality before the first observation was therefore unaccounted for, leading to some degree of overestimation. It was possible to estimate historical duckling survival in the Netherlands using the data provided by Eygenraam (1957), which amounted to a 43–53% survival rate from day one to 8 weeks for the years 1951–1954 ( Supplementary Material Table S2 (duac020_suppl_supplementary_material.docx)). Duckling survival thus appears to have halved since the 1950s. Similarly high values as reported by Eygenraam (1957) have been reported for Mallards in Canada (Bloom et al. 2012) and in South Dakota (Stafford and Pearse 2007), where adult survival and nest success appeared to be the factors limiting population growth (Hoekman et al. 2002). Duckling survival rates similar to the current ones in the Netherlands were found in Southland, New Zealand (0.16–0.30; Sheppard 2017), the Prairie Pothole Region in North Dakota (0.18–0.22; Amundson and Arnold 2010, 2011), and Central Valley, California (0.25; Chouinard Jr and Arnold 2007). For these regions, rapid population declines were observed (λ = 0.74–0.84) and agree with our sensitivity analysis and LTRE in pointing towards duckling survival as the limiting factor for the observed combination of vital rates (Cowardin et al. 1985, Amundson et al. 2013, Sheppard 2017).

More evidence that duckling survival limits population growth in the Netherlands comes from the estimated duckling survival rate for the years 2003–2017, acquired through the IPM: a mean duckling survival rate as low as 18% was required to explain the population trend given the observed estimates of the other vital rates. Although this duckling survival estimate should be interpreted with caution, as there were no explicit data for these years, our findings do indicate that the other vital rates were together unable to explain the population decline observed.

FIGURE 3.

Annual estimates of clutch size (A), nest success (B), duckling survival (C) and adult survival (D) for Mallards in the Netherlands based on an integrated population model. Error bars around annual estimates represent 95% Bayesian credible intervals (CI). Solid and dotted horizontal lines depict the parameter mean across all years and the corresponding CI.

TABLE 1.

Values used in, and the results of, a life-stage sensitivity analysis and a life-table response experiment (LTRE) for the Dutch population of the Mallard for the years 2003–2020. Historical sources: (1) Clutch size and duckling survival for 1951–1954 estimated (see  Supplementary Material Table S2 (duac020_suppl_supplementary_material.docx)) using data from Eygenraam (1957); (2) Adult survival for 1911–1965 as reported by van den Bremer et al. (2015).

Clutch size. Clutch size was stable throughout the study period and was not correlated with estimates of λ. This is in agreement with sensitivity analyses which suggest that clutch size is a rather insignificant predictor of population growth (Hoekman et al. 2002, Sheppard 2017). Clutch size estimates were surprisingly low across all years with on average 8.2 eggs per nest. Eygenraam (1957) reported an average brood size for one-day-old ducklings of 9.2–10.7 for 3,717 broods in the period 1951–1954 in the Netherlands, indicating a mean clutch size of at least that many eggs. Reported clutch sizes from the USA, Canada and New Zealand were all higher as well, with average clutches ranging from 8.3 up to 11.5 (e.g., Hoekman et al. 2002, Zicus et al. 2003, Devries et al. 2008, Sheppard 2017). In our study, clutch size might have been underestimated due to observed nests being predominantly from cultivated landscapes. Egg production in Mallards is known to be dependent on the nutrition available to the female (Krapu 1981, Pehrsson 1991) which varies between habitats (Ball et al. 2002). Agricultural grasslands may be of less nutritional value compared to natural systems. In the Netherlands, the availability and nutritional quality of important food items for many meadow birds are known to have decreased, due to different aspects of agricultural intensification, such as artificial lowering of water tables, fertilization, and increased mowing frequencies (Sanders et al. 2004). It should be noted that if clutch size has truly been underestimated in our study, the predicted duckling survival for 2003–2017 would in fact be even lower than reported here.

Nest success. Nest success throughout the study period was high and also the most variable of all low-level vital rates. Mean nest success in the Netherlands was considerably higher than estimates reported by large-scale North American studies, which often range between 10 and 25% (e.g., Klett et al. 1988, Davis 2008). Similar studies for European populations are lacking. High nest success in the Netherlands would explain the lack of a relationship between nest success and reproductive rate in this study, as other vital rates become much more limiting. For example, Hoekman et al. (2002) estimated that for Mallards, nest success only starts being a limiting factor for population growth at values below 20%.

Adult survival. We estimated a singular survival rate of 71% for all Mallards due to limited information on sex and age for many individuals. Previously reported sex and age effects in Europe have been of relatively minor importance, with differences in annual survival rate of 1-2% (Gunnarsson et al. 2008, 2012). We therefore do not expect our singular estimate to deviate strongly from reality. Although λ was sensitive to fluctuations in adult survival, this vital rate was stable throughout the study period and also higher than historical estimates for the Netherlands (67%; van den Bremer et al. 2015) and similar to recent estimates for northwest Europe in general (69–71%; Gunnarsson et al. 2012). Mean adult survival was also high compared to other countries (e.g., Hoekman et al. 2002, McDougall and Amundson 2017), possibly a result of relatively low hunting pressure in the Netherlands. Estimates for both nest success and adult survival, which were found to be the most important limiting factors to Mallard population growth in other studies, are therefore high and unlikely limiting to λ in the Dutch population (Hoekman et al. 2002). These findings, combined with low estimates for duckling survival, give a strong indication that the latter is the limiting vital rate to Mallard population growth in the Netherlands.

Management and Research Implications

The IPM revealed that duckling survival is likely the limiting factor to Mallard population growth, so management action to halt the population decline should primarily aim to increase this vital rate. However, why exactly duckling survival rates are so low and appear to have declined since the 1950s is not yet clear, making it difficult to design targeted measures. We identify two candidate drivers that may be able to explain the low duckling survival. First, food availability for ducklings may have been reduced due to improved water quality, which is likely to have reduced the bulk abundance of aquatic invertebrates benefitting from eutrophication (Schekkerman et al. 2016; Hallmann and Jongejans 2021). Secondly, the decreased duckling survival may be due to increased predation pressure. The populations of several important predators, most notably Lesser Black-backed Gull (Larus fuscus) and Common Buzzard (Buteo buteo), have strongly increased since the 1990s and might have increased predation of Mallard ducklings (Sovon 2019). Interestingly, the Dutch population of the Gadwall (Mareca strepera), an ecologically related dabbling duck, has increased eightfold since the 1990s (Sovon 2018). Competition between the two species is unlikely to have caused a decrease in Mallard duckling survival, since Mallards still outnumber Gadwalls by a factor of 10 (Sovon 2019). Although comparative studies between the two species are lacking in Europe, nest record data and field observations indicate that Gadwalls breed later and that Gadwall females with broods are more cryptic than Mallards (Schekkerman et al. 2016). Gadwall brood survival might therefore be higher due to relatively low predation compared to Mallard broods, which are not as well hidden by the female. Future studies should therefore focus not just on Mallard broods but on Gadwall broods as well.

This study shows that integrated population modeling can provide useful insight into population dynamics even if explicit data are partly missing. By incorporating the pieces of information that are at hand, the holistic and flexible nature of IPMs makes it possible to make rough estimations of missing vital rates. Nevertheless, the provision of real data remains crucial to the monitoring of waterfowl populations (Elmberg et al. 2006). The Mallard could be considered the best-studied duck species in the Netherlands, but even then data for migration, hunting bag, and reproductive rates outside farmland remain scarce. We therefore reiterate the plea by Elmberg et al. (2006) for a coordinated international effort to monitor demographic processes in dabbling ducks.