Large-scale changes in marine and terrestrial environments drive the population dynamics of long-tailed ducks breeding in Siberia

Migratory animals experience very different environmental conditions at different times of the year, i.e., at the breeding grounds, during migration, and in winter. The long-tailed duck Clangula hyemalis breeds in the Arctic regions of the northern hemisphere and migrates to temperate climate zones, where it winters in marine environments. The breeding success of the long-tailed duck is affected by the abundances of predators and their main prey species, lemmings Lemmus sibiricus and Dicrostonyx torquatus, whose population fluctuation is subject to climate change. In the winter quarters, long-tailed ducks mainly eat the blue mussel Mytilus edulis. We examined how North-west Siberian lemming dynamics, assumed as a proxy for predation pressure, affect long-tailed duck breeding success and how nutrient availability in the Baltic Sea influences long-tailed duck population size via mussel biomass and quality. Evidence suggests that the long-tailed duck population dynamics was predator-driven on the breeding grounds and resource-driven on the wintering grounds. Nutrients from fertilizer runoff from farmland stimulate mussel stocks and quality, supporting high long-tailed duck population sizes. The applied hierarchical analysis combining several trophic levels can be used for evaluating large-scale environmental factors that affect the population dynamics and abundance of migrants from one environment to another.


Materials and methods
Ethics declaration. This article does not contain any experiments on animal subjects performed by any of the authors.
Observations of long-tailed ducks. Spring migration counts of long-tailed ducks have taken place at Söderskär (60°07′ N, 25°24′ E, Fig. 1), Gulf of Finland each spring 1968-2014. These counts are generally assumed to reflect population size 8 , and there is no indication that the main migration corridor along the northern coast of the Gulf of Finland (Fig. 1) for Baltic long-tailed ducks has changed during the study. However, unfavourable wind direction and velocity can force migrating birds off their preferred migration route in the vicinity of the Söderskär observatory 52,53 , thereby causing variation in migrant bird counts. To control for this error source, we used daily wind direction and velocity measures at Söderskär during spring migration surveys. We developed a weighted wind direction and velocity factor to be implemented into the observation model of the long-tailed duck to correct long-tailed duck population size estimates by reducing noise due to the effect of wind on spring migration counts 54,55 (details in Supplementary Methods).
We estimated duck recruitment from the number of juvenile (1st to 2nd calendar year individuals) in relation to the number of adult female birds shot by hunters in autumn and winter (juvenile proportion in year t refers to bags from October in year t until February in year t + 1) from the Danish Wing Survey data 56,57 . We then related this to lemming abundance (of the year t-1) in North-western Siberia during the breeding period and climate conditions in breeding and wintering grounds (Supplementary Methods). Hunted long-tailed ducks in Danish waters belong to the North-west Siberian breeding population with the majority of long-tailed ducks wintering in the southern Baltic Sea 6,10 , and these birds migrate along the Gulf of Finland also including the Söderskär observatory 8,52 . Thus, we assumed that the Danish wing samples represent environmental impacts at the breeding grounds. These processes should be reflected to spring migration numbers as well.

Environmental variables.
(1) Nutrient levels. Dissolved inorganic nitrogen (DIN), considering the sum of the oxidized nitrogen and ammonium pool, and dissolved inorganic phosphorus (DIP) were estimated for the Baltic Sea major basins each year from 1970 to 2016 47 . As the variables for nutrient availability at the main wintering area of the long-tailed duck, we used the annual totals of DIN and DIP recorded from two major basins in the southern Baltic Sea, the Baltic Proper and the Danish Straits 47 (Fig. 1). Nutrients can have positive or negative effects on the long-tailed duck population size. Additional DIN stimulates primary production and thereby growth of mussels, the primary food of this species during winter [41][42][43] . Negative effects can arise because of hypoxia and bottom death that can occur at low DIN:DIP ratios 47 in areas of poor mixing 46 . The total annual amount of fertilizer applied by Danish farmers is a reliable proxy for total nitrogen and phosphorus runoff into the marine environment 48 60,61 (data downloaded from www. world clim. org) for the three regions in Western Taimyr for which lemming abundance data were available. Climate variation was averaged for the three areas corresponding to the three rodent surveys (Fig. 1). We estimated climate variability impacts on rodent population sizes 38 , which in turn may influence the breeding success and population dynamics of the long-tailed duck 8 . Similarly, for estimating climate-related consequences on long-tailed duck juvenile proportions, we averaged annual precipitation and temperature scores for North-western Siberia (Fig. 1). For processing the climatological map data, we used the R package 'raster' 62 . (4) North Atlantic Oscillation Index (NAOI 63 ). Winter NAOI (December in year t-1 until February in year t) is an index of the severity of winter conditions in the Northern Atlantic 63,64 . High index values indicate above average winter temperatures and high levels of precipitation in the Baltic Sea 64 . We investigated the association between winter NAOI and the proportion of juveniles killed by Danish hunters and ultimately the population dynamics of the long-tailed duck. (5) Abundance and quality of blue mussels. At the wintering grounds, long-tailed ducks feed on blue mussels, but no long-term data exist for mussel stocks in the Baltic Sea. Comparable data are available from the Wadden Sea, where annual blue mussel stocks were estimated for the intertidal zone of the Danish (1986-2007, and 2017) and Schleswig-Holstein (Germany, 1998-2015) parts of the Wadden Sea ( Fig. 1) during autumn by ground sampling on mussel beds and estimation of mussel beds from aerial photography 65,66 . Mussel stock data are available for both areas from 1998 to 2007 and annual biomass showed a positive correlation between the areas (r = 0.920, N = 10). Blue mussel quality data, measured as flesh to shell ratio 67,68 , were available for each autumn from 1998 to 2013 from 19 sites in the Baltic and Wadden Sea (see Supplementary Methods).
The region of the mussel surveys are outside the core winter grounds of the long-tailed duck 6 . Thus, we did not accept the mussel data in the ultimate model for long-tailed duck population dynamics. Instead, we aimed at indicating the relation between nutrient load (in the form of fertilizer use in the Danish farmland 58 ) and blue mussel stocks on the Wadden Sea 69 . A relation found between fertilizers and mussels on the Wadden Sea would indicate that the availability of dissolved nutrients (DIN and DIP as linked to fertilizer use 58 ) can be used as a proxy for mussel abundance and quality in the southern Baltic Sea. The results shown for other seabird populations in that region 70,71 support this reasoning.
We also aimed at showing how winter and spring temperatures affect blue mussel stocks on the Wadden Sea. For these, we calculated regional temperatures using the same global climate database as described above (see "North-west Siberian climate"). In the Baltic Sea, in the western part of the Finnish Archipelago, blue mussels grow from larvae to 10 mm within 3-4 years 72 so we expect mussels in the southern Baltic Sea ecosystem to attain the size of 9 mm preferred by long-tailed ducks 44 in about the same period after spawning. Hence, assuming that increasing winter temperatures reduce mussel quality 68 and energy for spawning 70 , a reduction in 9-mm mussel cohorts should emerge a few years after mild winters.
Hierarchical modelling. We used an integrated hierarchical model 73 , which is a complex stochastic system partitioned into a dependent sequential set of simpler sub-models dynamically affecting the performance of the main system of interest 73 (Supplementary Methods, Supplementary Table S1).
Lemming population dynamics z (t,j) , combining three separate population time-series j of lemming abundances in the Western Taimyr Peninsula (1965-2017, t = 0 refers to 1965), were estimated based on the following state-space model: in which random variation in lemming population size estimates, z (t,j) , are negative binomially distributed around the expectation, z (t,1) (Eq. (1), detailed in Supplementary Methods). a z is the intercept for j = 1 , and c z is a density dependent parameter 51 . The parameter b (C)k controls for the response of population j = 1 to climate variation ( C (t,k,j=1) ) . Subscript k indicates precipitation in summer (t, k = 1) and the previous autumn (t − 1, k = 2) for each population j . Parameter vectors a and b measure the consistency in three-year cyclic dynamics expressed in dummy variable (0, 1) matrices X 1 (1965-1994) and X 2 , (1995-2017) in which the number one enables regular effects by three-year intervals (details in Supplementary Methods). b t is a parameter for a temporal trend. ǫ z(t,j) controls for normally distributed environmental disturbances. Additive (to z (t,1) ) scaling parameters for populations 2 and 3 are denoted by δ 2 and δ 3 ; similarly, climate effects for populations 2 and 3 are quantified with θ k,2 and θ k, 3 .
We constructed a logit model, Eq. (2), to generalize the variation of juvenile proportions during the period 1967-2017, based on the number of juveniles per adult female shot by Danish hunters: (1) www.nature.com/scientificreports/ Here, p (t) is the expected proportion of juveniles (in autumn and winter). a p is the intercept. The parameter b L controls for the responses of p (t) to estimated (log) lemming population size in the previous year based on the longest series of lemmings (Kara Sea, z (t−1,j=1) ) from Eq. (1). b C parameterizes the variation related to winter climate C (t−3) (NAOI Dec-Feb ) three years before. The terms b (S,p)k S (t,k) quantify the effects of precipitation (k = 1) and temperature (k = 2) on the expected juvenile proportions in North-western Siberia in May and June, and these were weighted, e.g. for precipitation (P), as S t,1 = 1 − w June,1 P May + w June,1 P June . The weight parameter, w June,1 , is beta-distributed varying between zero and one (uninformative prior beta-distribution, Beta(1, 1) 74 ; Supplementary Methods) r p(t) controls for random effects (full description in Supplementary Methods).
A log scale state-space model, Eq. (3a), quantifies fertilizer effect on nutrient pools measured in the Danish Straits j = 1 and the Baltic Proper j = 2 , returning state estimates D (t,k) for DIN (k = 1) and DIP (k = 2) for the southern Baltic Sea in 1970-2016: where a k parameters are intercepts for DIN and DIP, and c k controls for serial correlation. b k is a parameter for fertilizer-use (F t ) effect on the nutrient pools. Non-linear effects are specified with a spline function for fertilizers generating two additive smoothing variables 75 : s(F) t,k . The lag parameter is binomially distributed and can have values of 0, 1, 2 or 3 (years) (Supplementary Methods). e (t,k) is a normally distributed random effect term (Eq. (3a)). Observed nutrient amounts link to the states via an observation model, Eq. (3b), with intercepts a j,k and standard deviations s D(j,k) . In the following system, the state-space model, Eq. (4), expressing log-scale long-tailed duck dynamics n (t) (1968-2014) is written briefly as: in which n (t) is the expected population size in year t . n (t−1) is the population size estimate of long-tailed ducks in year t − 1 that has a negative binomial error structure 50,76 around the expected population size n (t−1) . a is the intercept, and c is a density-dependent parameter whose value can vary from zero to one. R (t−1) = 1 + p (t−1) , adjusted by the parameter β R , is a variable for recruitment based on estimated juvenile proportion ( p (t−1) ) in autumn and winter preceding spring migration from Eq. (2). Dissolved nitrogen DIN (t−1) and phosphorus DIP (t−1) describe annual nutrient pools in the southern Baltic Sea (outputs from Eq. (3a)). The nutrient series implemented in Eq. (4) were scaled to have a mean of zero and variance of one and modelled as fixed environmental variables. w D weights nutrient effects. Parameters β R and w D are beta-distributed varying between zero and one (uninformative prior beta-distribution, Beta(1, 1) 74 ). The parameter β D is normally distributed. ε (t) controls for random environmental disturbances (Supplementary Methods). The observation model for Eq. (4) is specified with Gaussian errors and a priori assumed environmental disturbances (wind) as: y (t) |n (t) ∼ normal n (t) + γ x (t) , τ 2 , where γ is a parameter for east-west aspect winds x (t) , and τ is the standard deviation of a random observation-error process (Supplementary Methods). The joint-likelihood 73 of the integrated hierarchical model combining Eqs. (1), (2) and (4) is summarized with Supplementary Table S1.
To estimate how biomass of mussels increases with increasing nutrient availability, we built a model based on two mussel biomass surveys of the Wadden Sea, assuming that temperature and nutrients have comparable effects on mussel population growth in the Wadden Sea and the southern Baltic Sea (Supplementary Methods). These effects exclude hypoxia due to nutrient overflow because this does not occur in the strongly tidal wellmixed Wadden Sea, though it does in the Baltic Sea. Thus, a state-space model for mussel (log) biomasses in The Danish Wadden Sea ( j = 1 ) for the period 1986-1997 is expressed as: where a µ is the intercept and c µ is a density-dependent parameter with the term b F F t−lag µ controlling for fertilizer effect on biomass estimates µ (t,j=1) . s(F) t−lag µ is an additive smoothing variable 75 as explained for Eq. (3a). The parameter lag µ can have values 1, 2 or 3 (cf. Eq. (3a)). The function f (T) parameterizes the weighted effects of mean winter temperature (December-February) on biomass estimates two and three years later. e µ(t) controls for state-process errors due to environmental variance and δ µ(t,1) controls for demographic variance 77,78 . Biomass estimates µ (t,j=2) for Schleswig-Holstein follows the parameterization given in Eq. (5), excepting an independent demographic error term δ µ(t,2) and additive scaling parameter for an observation model (Supplementary Methods). A logit model for flesh proportion in mussels f (t) for the period 1998-2013 was expressed as:

Results
The magnitude of lemming population peaks in the western Taimyr Peninsula have declined over time (Fig. 2) and the 3-year cyclic regularity was broken after 1994 (Supplementary Colder and wetter springs and summers in North-western Siberia were followed by fewer long-tailed duck recruits on the wintering grounds in the Baltic. The precipitation variable (Fig. 3a), with weighting particularly   Table S1b), had a stronger effect than the temperature variable (Fig. 3b).
Together the North-west Siberian climate variables explained 26% of the total variation in juvenile proportions. The proportion of juveniles declined the year following lemming peaks (Fig. 3c), implying that lemming population troughs that followed high abundance peaks were poor years for long-tailed duck recruitment. Lemming population variation explained 9% of juvenile proportions (Supplementary Table S1b). High NAOI was followed by declines in juvenile proportion with a three-year lag (Fig. 3d), suggesting poor recruitment three years after a warm winter. Together, the climate parameters (Fig. 3) explained 39% of the total variance estimated for the long-tailed duck juvenile proportion (Supplementary Table S1b). Long-tailed duck spring migration counts as estimated with Eq. (4) increased with increasing DIN and decreasing DIP in the southern Baltic Sea in the previous year (24% variance partition proportion for the joint effect with high statistical confidence, Supplementary Table S1c). The weighting parameter w D (= 0.49) suggested similar effect sizes for DIN and DIP. Figure 4a shows the joint effect. Increases in juvenile proportions in autumn and winter, i.e. R (t) from Eq. (2), were followed by slight increases in numbers of spring-migrating long-tailed ducks n (t+1) from Eq. (4) (5% variance partition proportion, Supplementary Table S1c). This trend is best illustrated by plotting annual population growth rates, n (t+1) − n (t) , against R (t) , although this relationship is noisy (Fig. 4b,c) with parameter distribution skewed rightward: P(β R > 0.05) = 0.9 (mean of β R = 0.25 ± 0.17 SD). Density-dependent population regulation of the long-tailed duck was observed (22% variance partition proportion).
The observation process variation (i.e., sampling error or observation error) of spring migration counts was affected by east-west aspect winds such that observed counts decreased with strengthening west winds during the annual census period. The details of wind effects on the observations, which delineated the performance of the population model given in Eq. (4), are reported in Supplementary Methods. After correcting for measurement www.nature.com/scientificreports/ errors, the number of spring migrant long-tailed ducks showed relatively smooth annual variation and increased from the end of the 1960s up to 1991 followed by a steep decline from 1992 until the end of the 2000s. Some recovery took place in the beginning of the 2010s (Fig. 4d). The whole system with respect to the long-tailed duck is summarized in Fig. 5. Application of more fertilizer to Danish farmland led to higher DIN and DIP in the southern Baltic Sea the following year. This 1-year lagged response (lag = 1) resulted from 94.25% of the simulation updates (Supplementary Table S2).
High winter temperature was followed by lower mussel stocks in the Wadden Sea 2-3 years later (more weight for 2-year lag). Mussel biomass increased with the increasing use of fertilizers (Fig. 6a) in previous years, most strongly 1-2 years before. Simulation chains revealed frequencies for the lags of 0, 1, 2, and 3 years of 6.3, 31.0, 43.4, and 19.3 percent, respectively. Fertilizer and the smoother based on the spline function dominated the estimated variance for mussel biomass (Supplementary Table S3a). More demographic variance caused by local disturbances affecting mussel dynamics was found for the population in the Danish than the Schleswig-Holstein part of the Wadden Sea (Supplementary Table S3a). Population trends were not perfectly parallel in the two areas (Fig. 6b). Together, mussel biomass showed large variation during the period 1986-1998, entering steep Lemmus sibiricus and Dicrostonyx torquatus abundances and juvenile proportions of LtD. "Precipitation" and "Temperature" refer to North-western Siberia in June and July, and "NAOI" is the North Atlantic Oscillation Index. Precipitation in the Western Taimyr Peninsula affects lemmings. Red/blue arrows indicate a negative/ positive effect. Solid lines indicate direct effects. Indirect effects (long dash lines) involve trophic cascade forcing from one trophic level to another. Mechanisms behind Siberian climate effects (short dash lines) on lemmings and juvenile proportions were out of the scope of this study and were not included in the general discussion. Black arrows denote random errors from unknown sources and circle arrow illustrates density dependence. Percentages for DIN and DIP describe direct fertilizer plus smoother effects based on the fertilizer series. Variance partition proportions (%) shown for the predictors. The sampling error proportion of the total LtD system (including sampling error variance) is in parentheses. Long-tailed ducks (female and hatchlings, and male and female on the sea), lemming L. sibiricus, and fertilizer-use symbol (crop) drawings by Kati Rintala. www.nature.com/scientificreports/ declines from the turn of the 1990s and recovering from 2011 up to 2017 (Fig. 6b). The quality of blue mussels in autumn in terms of (log) flesh/shell ratio of the total biomass, i.e. logit f p(t) from Eq. (6), increased with the previous winter's temperature (19% variance partition proportion, Supplementary Table S3b) and fertiliser use over the two previous years (42% variance partition proportion, Supplementary Table S3b, Fig. 7a,b) but spring temperatures had no effect on mussel quality. These results suggest that mussels in the Wadden Sea were in better www.nature.com/scientificreports/ condition, as measured in autumn, during the years of high nutrient availability and after warm winters. The system with respect to blue mussel flesh/shell ratio and biomass is summarized in Fig. 8.

Discussion
This integrated hierarchical model for the long-tailed duck population investigates a complex stochastic system affected dynamically by sub-systems that are also stochastic and regulated by independent variables 73,79 . This model allowed us first, to correct for population estimate errors due to variation in wind velocity and direction at Söderskär, Gulf of Finland, during spring migration, which accounted for 51% of the total observation error variance. This result underlines that ignoring wind effects on counts would have led to serious misspecifications in the state-space model for long-tailed ducks. This model allowed us to identify and quantify the variables contributing to population size variation of the long-tailed duck by using sub-models for trophic processes driven by nutrient and climate factors occurring at both the breeding and the wintering grounds. The model attributes 51% of total variation in long-tailed duck population size, estimated from spring migration counts, to ecological processes occurring at the breeding or overwintering grounds (Fig. 5). The processes influencing long-tailed duck populations in these two habitats hypothetically are driven by different ecological interactions, for example predator-driven on the breeding grounds 12,14,80 and resource-driven on the wintering grounds 45,58,70,71 . For example, precipitation evidently influences lemming dynamics on the breeding grounds while temperature and nutrient runoff from the land are thought to affect mussel reproduction, growth and survival on the wintering grounds.
On the breeding grounds, lemming population cycles were related to the number of new recruits of long-tailed ducks, estimated from the proportion of juvenile birds killed by hunters in the winter. The year following high lemming years had few recruits while those following low lemming years had many recruits (Fig. 3c). This is as expected if predator numbers are high the year after a lemming peak and if these predators switch to eggs and nestlings when lemmings are rare 12,16,81 . Western Taimyr lemming dynamics was associated with, and presumably driven by climate, particularly autumn and winter precipitation 38 (Fig. 5). High autumn precipitation was associated with low lemming population size the next year. Heavy autumn precipitation may cause ice formation at the soil surface in the subnivean space, reducing the insulation properties of snow 39,40 , which leads to poor lemming winter success 22,38,[82][83][84] .
High late spring and early summer precipitation with low temperatures in North-western Siberia were related to decreases in juvenile proportions in long-tailed ducks (Fig. 5). This may be linked to varying nesting areas and predation pressure during snow melt: more precipitation during cold springs delays snow melt 85,86 , which limits open areas for nesting and thus enhances predation on birds' nests during early breeding season 86 .
Recruitment, as estimated from the juvenile proportion of hunters' returns, explained 5% of variation in spring migration numbers. This weak association is as expected for a long-lived species like the long-tailed duck 87 , where adult survival drives population dynamics more than does fecundity [88][89][90][91] . Furthermore, the juvenile proportion of hunters' returns is probably representative of the population. Long-tailed duck juveniles resemble www.nature.com/scientificreports/ mature females in autumn and winter so hunting pressure is unlikely to be age dependent, as for some other quarry species 92 .
On the wintering grounds long-tailed ducks feed mainly on blue mussels 43,44,68 . Mussels reproduce more during cold than warm winters 45,68 . NAOI values that represent variation in winter temperatures showed the strongest association with juvenile long-tailed duck proportion at the wintering grounds in the southern Baltic Sea three years later. If mussels require two or three years to reach the preferred size for long-tailed duck food items 44,72 , then food on the wintering grounds will be scarce two or three years after a warm winter but abundant two or three years after a cold winter. In the Wadden Sea the best-fit lag between cold winters and mussel biomass increase was two years (Fig. 8). Years of high food abundance on the wintering grounds should enhance female condition, increasing their fecundity and maternal investment 45,70 . Thus, we would expect, as we find, higher nesting and fledging success three years after cold winters.
We detected apparent consequences of fertilizer runoff on long-tailed duck populations via the bottomup forcing primed by dissolved nutrients 47 . Long-tailed duck populations increased with increasing DIN that, through increasing phytoplankton abundance, likely translated into increased biomass and quality of mussels. However, the opposite trend was observed for DIP because increased DIP led to larger and more longlasting hypoxia and bottom death 46,47,93 . Similar results have been seen in related studies of eiders and other waterbirds 45,58 .
Here we estimated the bottom-up forcing of dissolved nutrients, in the form of nutrient load from Danish farmland, on mussel biomass data from the Wadden Sea. Though we have argued that the predictions about how DIN influences food supply for long-tailed ducks are likely to be appropriate, the results on the Wadden Sea mussels related to DIP cannot be directly generalized to the Baltic Sea. This is because, unlike the Wadden Sea, the Baltic Sea experiences hypoxia due to nutrient overload affecting the cycling of phosphorus 46,47,93 . Hypoxia in the Baltic Sea has resulted in the elimination of benthic fauna over large areas and has severely disrupted food webs 58,94-97 with expectable adverse effects on sea birds 6,58,98,99 . It is known that hypoxia increases with increasing DIP in the Baltic Sea 47,93 , which accords with the apparent decline in long-tailed ducks abundance with increasing DIP.
In conclusion, we used hierarchical Bayesian models to quantify the effects of a complex array of processes on the abundance of the long-tailed duck, an Arctic migratory bird. This approach allowed us to estimate the amount of variance accounted for by environmental conditions on both the breeding and wintering grounds. Considering drivers at the Arctic breeding grounds, our analysis emphasizes the hypothesis that predation limits nesting success when lemming predators shift to alternative prey such as the long-tailed duck. In the marine winter quarters fertilizer runoff from agriculture apparently drives a bottom-up effect, stimulating food availability by increasing primary production and thereby mussel biomass and quality. However, phosphorus overload could limit food availability by inducing hypoxia and bottom death. These conditions at the wintering grounds affect juvenile survival and adult female body condition before spring migration to the Arctic breeding grounds. The originality of this study resides in our ability to quantify the effects at both breeding and wintering areas that present very different ecological and environmental conditions and challenges. This approach should be useful for analysing the dynamics of other migratory species confronted with divergent environmental conditions and trophic interactions in varying habitats.

Data availability
The data supporting the results of this study are provided as Supplementary Data (R workspace list-type objects printed as "data_*.docx" and "inits_*.docx" files) referenced in the Bayesian analyses (section "Program codes" in Supplementary Methods).