Integrated population modeling identifies low duckling survival as a key driver of decline in a European population of the Mallard

ABSTRACT Europe's highest densities of breeding Mallards (Anas platyrhynchos) are found in the Netherlands, but the breeding population there has declined by ∼30% since the 1990s. The exact cause of this decline has remained unclear. Here, we used an integrated population model to jointly analyze Mallard population survey, nest survey, duckling survival, and band-recovery data. We used this approach to holistically estimate all relevant vital rates, including duckling survival rates for years for which no explicit data were available. Mean vital rate estimates were high for nest success (0.38 ± 0.01) and egg hatch rate (0.96 ± 0.001), but relatively low for clutch size (8.2 ± 0.05) compared to populations in other regions. Estimates for duckling survival rate for the three years for which explicit data were available were low (0.16–0.27) compared to historical observations, but were comparable to rates reported for other regions with declining populations. Finally, the mean survival rate was low for ducklings (0.18 ± 0.02), but high and stable for adults (0.71 ± 0.03). Population growth rate was only affected by variation in duckling survival, but since this is a predominantly latent state variable, this result should be interpreted with caution. However, it does strongly indicate that none of the other vital rates, all of which were supported by data, was able to sufficiently explain the population decline. Together with a comparison with historic vital rates, these findings point to a reduced duckling survival rate as the likely cause of the decline. Candidate drivers of reduced duckling survival are increased predation pressure and reduced food availability, but this requires future study. Integrated population modeling can provide valuable insights into population dynamics even when empirical data for a key parameter are partly missing. LAY SUMMARY Demographic data and population modelling are pivotal for identifying limiting processes and vulnerable life stages involved in the decline of wildlife populations. The Mallard (Anas platyrhynchos) population in the Netherlands has decreased by 30% since 1990 but the causes remain unclear. We combined the available historical and current data for Mallard population size, reproduction, and survival of ducklings and adults in a population model to investigate the processes involved in the decline. Duckling survival was the only vital rate to strongly affect variation in population growth and was low compared to both historical data for the Netherlands and to duckling survival estimates for stable Mallard populations in other countries. Future research should investigate why duckling survival is low and how conservation action can target this life stage. SAMENVATTING In Nederland vindt men de hoogste dichtheden broedende Anas platyrhynchos van Europa, maar sinds de jaren 90 is deze populatie met 30% afgenomen. De oorzaak hiervoor is echter onduidelijk. In dit onderzoek worden recente data over de populatieaantallen, nestgegevens, kuikenoverleving, en adulte overleving van Anas platyrhynchos gecombineerd in een geïntegreerd populatiemodel. Deze allesomvattende methode is gebruikt zodat alle relevante demografische processen vastgesteld kunnen worden, waaronder de kuikenoverleving in jaren waar geen expliciete data voor beschikbaar waren. De jaarlijkse overleving van adulten bleek stabiel en hoog (0.71 ± 0.03) tijdens de studieperiode, terwijl het nestsucces hoog (0.38 ± 0.01) maar variabel was. De gemiddelde legselgrootte was relatief laag (8.2 ± 0.05). De kuikenoverleving voor de drie jaar waar data voor aanwezig waren bleek laag (0.16–0.27) in verhouding tot historische waardes, maar was vergelijkbaar met waardes zoals geobserveerd in andere regio's met afnemende populaties. Voor de hele periode werd de jaarlijkse kuikenoverleving op 18% geschat. Kuikenoverleving was de enige parameter die correleerde met de populatiegroeisnelheid, hetgeen echter voorzichtig geïnterpreteerd moet worden omdat kuikenoverleving grotendeels een latente variabele was. In combinatie met een vergelijking met historische waardes wijzen onze bevindingen echter toch op een afgenomen kuikenoverleving als een waarschijnlijke reden voor de populatieafname van Anas platyrhynchos. Mogelijke oorzaken hiervoor zijn toegenomen predatiedruk en afgenomen voedselbeschikbaarheid, maar hier is vervolgonderzoek voor nodig. Dit onderzoek laat zien hoe geïntegreerde populatiemodellen waardevol inzicht kunnen bieden in de populatiedynamiek van een soort, zelfs als empirische data voor bepaalde jaren ontbreekt.


INTRODUCTION
A solid understanding of the interplay between individual performance and population growth is crucial in order to implement effective measures to halt wildlife population declines. Population modeling can provide estimates of these demographic processes and reveal which ones are limiting population growth, indicating which life-history components should be targeted by conservation action (Mills 2012). Ideally, this requires long-term, standardized monitoring schemes to directly track changes in both population abundance and vital rates over time (Elmberg et al. 2006). For waterfowl, this includes at least data for population size, adult survival, and reproductive rates, along the entire range of countries that are part of a species' flyway (Elmberg et al. 2006). Although monitoring programs are increasingly put into place (e.g., Du Feu et al. 2016), a complete picture of demographic rates remains rare, sometimes even for abundant species. Historical records of baseline vital rates are often unavailable as well, making it difficult to directly attribute population declines to a decrease of a specific vital rate. Despite its abundance, the Mallard (Anas platyrhynchos) in the Netherlands is no exception. It is the most numerous breeding waterbird species and the second most common non-passerine in the country, but the currently available demographic information has proven insufficient to explain the population decline in recent decades (Boele et al. 2014, van den Bremer et al. 2015. The Mallard is likely the most successful duck species in the world in terms of abundance and distribution (Keller et al. 2020). It plays an important role in seed dispersal between wetlands, and with an annual hunting bag of 4.5 million individuals, it is also the most hunted waterfowl species on the continent (Hirschfeld and Heyd 2005, Gunnarsson et al. 2012, Kleyheeg et al. 2017). In the Netherlands, the breeding population numbers around 200,000-300,000 breeding pairs and the winter population around 600,000-800,000 individuals (Sovon 2018). The Netherlands, therefore, houses the highest breeding densities of Mallards in Europe (Keller et al. 2020) and a significant proportion of the northwest-European flyway population, which numbers 4.2-6.7 million individuals (Nagy and Langendoen 2020). The Mallard breeding population in the Netherlands has been declining steadily across all habitats since national-scale monitoring of the population started in the 1980s. The decline has accelerated in the last decade, corresponding to a 30% population decrease between 1990 and 2019 (van den Bremer et al. 2015, Schekkerman et al. 2016, Sovon Vogelonderzoek Nederland 2019. A similar decline has not been observed in neighboring countries or in Europe in general, with a slight population increase since 1980 (Keller et al. 2020).
Population decline can be attributed to either a reduction in per capita survival and/or in reproductive rate, or to a net emigration of individuals outside of the area of interest. Long-distance migration data for the population in the Netherlands are scarce, but emigration is expected to play a minor role in Mallard population dynamics as females are highly philopatric (Oring and Sayler 1992, Blums et al. 1996. Long-distance breeding dispersal occurs only sporadically, and birds banded in the breeding season in the Netherlands are only rarely reported in another country during a following breeding season (Blums et al. 2003;https://vogeltrekatlas.nl/). The vast majority of Dutch Mallards are therefore likely year-round residents Ornithological Applications 124:1-12 © 2022 American Ornithological Society (Buijs and Thomson 2001). This is supported by a similar rate of decline of the winter population as observed in the breeding population (van den Bremer et al. 2015). The cause of the population decline is therefore most likely a domestic issue, and related to either reduced survival or reproductive rates, or both.
Sensitivity analyses of population growth rates in North American Mallard populations have revealed that population growth is particularly sensitive to 3 vital rates: nest success, adult survival, and duckling survival (e.g., Hoekman et al. 2002, Sheppard 2017. If emigration indeed does not play an important role, the most likely mechanism of population decline is a reduction in reproduction or survival. While data for adult survival and nest success are readily available, this is not the case for duckling survival, which is difficult to incorporate in standard nest monitoring schemes due to strongly precocial hatchlings. Since a study in the 1950s (Eygenraam 1957), data on Mallard duckling survival in the Netherlands had not been collected until a citizen science project to study ducklings was initiated in 2016, yielding sufficient data for analysis since 2018 (Kleyheeg 2019). No explicit data are therefore available for the main period of the population decline, making historical comparisons difficult.
The aim of this study was to estimate annual vital rates for the Dutch Mallard population for the period 2003 to 2020, in order to identify the demographic processes most likely responsible for the population decline. For this purpose, we constructed a Bayesian integrated population model (IPM) that combines multiple data sources, all gathered by citizen scientists, including those for vital rates and population abundance trends, to estimate the demographic processes of the Mallard since the turn of the century. The advantage of IPMs over more conventional population models is that datasets are not analyzed separately but are combined to form a single joint likelihood instead (Schaub and Abadi 2011, Arnold et al. 2018, Gamelon et al. 2021; Figure 1). Using this approach, demographic parameters can be estimated more precisely, and the uncertainty across all data sources is accounted for (Schaub and Abadi 2011). A third benefit is that the combined analysis makes it possible to estimate demographic vital rates for which empirical data are lacking (Schaub and Abadi 2011, Kéry and Schaub 2011, Nuijten et al. 2020. This feature of the IPM framework was pivotal for estimating duckling survival to fledging age for the years 2003-2017, which lacked empirical data on this specific life stage.

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.  (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 andThomson, 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(Gunnarsson et al. , 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 Ornithological Applications 124:1-12 © 2022 American Ornithological Society (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 1 ,..., y 20 and assume that each observation y t is made with constant variance σ y 2 . 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 ), 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 1 is equal to breeding incidence while p 2 through p 6 are described by p z = max{1 -0.26(z -1), 0}, where z is the z th nest attempt . 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 postfledging until the end of the breeding season was negligible. Observations of brood size (y b ) from a particular female at day τ 2 were modeled as a Poisson distribution with mean brood size reported during the previous observation of the female at day τ 1 . This was then multiplied by the daily duckling survival rate (φ d ) and raised to the power of the time interval between the two observations (τ 2 -τ 1 ): 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). 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 τ 1 and τ 2 , 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: Age=1 e βD0,t+βD,t * Age 1 + e βD0,t+βD,t * Age 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: 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 σ 2 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 timedependent 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). This yielded the following formulas: ϕ , and ε r,t = N 0, σ 2 r . 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 β D0,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 1 ∼ 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).  Table S1). Predicted values for the unobserved true population index N t conformed well with the observed annual index values (Supplementary Material Figure S1). 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). 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). 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).

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 λ.

DISCUSSION
To investigate the recent decline of the Mallard population in the Netherlands, we used an integrated population model to estimate all relevant demographic rates for the period 2003-2020 and their effect on the population growth rate. Estimates of clutch size and adult survival were stable, while those for nest success and duckling survival were not. The analyses revealed that λ was the most sensitive to changes in duckling survival, which was also the only vital rate strongly affecting reproductive rate and the vital rate with the largest LTRE contribution. These findings identify reduced duckling survival as a likely driver of the decline of the population since the 1990s.

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). 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 . 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; Arnold 2010, 2011), andCentral 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. 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 , 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 . 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(Gunnarsson et al. , 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 andAmundson 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 . 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 Ornithological Applications 124:1-12 © 2022 American Ornithological Society 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.

CONCLUSION
We estimated the relevant vital rates for the declining population of the Mallard in the Netherlands with an integrated population model using data from four different citizen science programs. Annual adult survival was high and stable while annual nest success was high and variable. Clutch size was relatively low compared to other reported estimates. Estimates of annual duckling survival could be calculated for 2018-2020 and were at levels low enough to limit population growth. Additionally, the model estimated duckling survival for the years 2003-2017 for which no empirical data were available using the holistic nature of integrated population modeling. Population growth was strongly correlated with reproductive rate, which in turn was correlated with duckling survival. This indicates that during the study period, no other vital rate than duckling survival was able to explain the population decline. Together with the comparison with historic values, these results are evidence that the Dutch Mallard population has declined since the 1980s due to a low duckling survival rate. This study shows that even in the face of missing data, integrated population modeling can provide useful insight into the factors that are limiting population growth. We recommend that the cause of low duckling survival rates in the Netherlands is investigated and emphasize the need for long-term, standardized monitoring projects of vital rates.