## 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 for Baltic long-tailed ducks has changed during the study. However, wind direction and velocity at Söderskär can force migrating birds off their preferred migration route 51,52, 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 counts. 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 53,54 (details in Supplementary Methods).

We estimated 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 55,56 and related this to lemming abundance (of the year *t*–1) in North-western Siberia during breeding and climate conditions in breeding and wintering grounds (Supplementary Methods).

## Environmental variables

(1) Nutrient levels. Dissolved inorganic nitrogen (DIN), considering the sum of the oxidized nitrogen and ammonium pool, and phosphorus (DIP) have been estimated for the Baltic Sea major basins each year from 1970 to 2016 46. We tested the effects of total DIN and DIP recorded from the southern Baltic Sea, the Baltic Proper and the Danish Straits major basins 46 (Fig. 1), covering the main wintering area of the long-tailed duck, on the population dynamics of this species. Nutrient effects can be positive, because additional DIN stimulates primary production and thereby growth of mussels, the primary food of this species during winter 40–42, or negative, because of hypoxia and bottom death that can occur at low DIN:DIP ratios 46 in areas of poor mixing 45. We used the total annual amount of fertilizer applied by Danish farmers during 1965–2016 as a predictor of marine nutrient pools in the southern Baltic Sea because this variable is a reliable proxy for total nitrogen runoff into the marine environment 47 (see Supplementary Methods).

(2) Lemming abundances. The majority of long-tailed ducks wintering in the Baltic Sea originate from the part of North-western Siberia including the Yamal and Taimyr Peninsulas 57. Lemming abundance across the entire breeding area should influence the dynamics of the long-tailed duck population, but such data are lacking. However, information on lemming abundance exists from three separate, long-termed surveys in the Western Taimyr Peninsula, North-western Siberia, which allow us to quantify regional population changes in lemmings (Fig. 1). One survey was performed within a zone of ca. 100 km along the coast of the Kara Sea, Western Taimyr, between 80 and 90°E 14,15 and two other surveys took place on the Western Taimyr Peninsula 13 (survey no. 44, Meduza Bay, 73°04´N, 80°30´E; and survey no. 45, Mys Vostochnyi, 74°06´N, 86°48´E) (see Supplementary Methods).

(3) North-west Siberian climate. We extracted annual (1965–2017) monthly climatological variables (precipitation and temperature) from a database of high spatial resolution global weather and climate data 58,59 (data downloaded from www.worldclim.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 tested if climate variability affected rodent population sizes 37, which may influence the breeding success and population dynamics of the long-tailed duck 8. Similarly, for estimating climate effects on recruitment, 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’ 60.

(4) North Atlantic Oscillation Index (NAOI). 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 61,62. High index values indicate high winter temperatures and high levels of precipitation in the Baltic Sea 61. NAOI data were obtained from https://climatedataguide.ucar.edu/climate-data/hurrell-north-atlantic-oscillation-nao-index-pc-based. We estimated the effects of winter NAOI on the proportion of juveniles killed by Danish hunters and on population dynamics of long-tailed ducks.

(5) Abundance and quality of blue mussels. Long-tailed ducks feed on blue mussels. Though no long-term data exist for mussel stocks in the Baltic Sea, such data are available from the Wadden Sea although this is outside the core winter grounds of the long-tailed duck 6. Two arguments justify the use of mussel stock data from the Wadden Sea as a proxy for the missing mussel stock data from the Baltic: First, population trends of mussel stocks are comparable over areas covering several hundred kilometres 63, and, second, the eider *Somateria mollissima*, another sea duck, also feeds primarily on blue mussels and winters in the Baltic Sea. Winter body condition of eiders in the Baltic is positively correlated with mussel stock estimates for the Wadden Sea 64,65.

In the Baltic Sea, in the western part of the Finnish Archipelago, blue mussels grow from larvae to 10 mm within 3–4 years 66 so we expect mussels in the southern Baltic Sea ecosystem to attain the size of 9 mm preferred by long-tailed ducks 43 in about the same period after spawning. 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 67,68. Mussel stock data are available for both areas from1998–2007 and annual biomass showed a positive correlation between the areas (*r* = 0.920, *N* = 10).

Data were available for blue mussel quality, measured as flesh to shell ratio 69, each autumn from 1998 to 2013 from 19 sites in the Baltic and Wadden Sea. The reduction in mussel flesh content during winter (October–March) due to temperature-dependent respiration was estimated for the Baltic Sea 70. This allowed us to estimate long-tailed duck food resources by predicting mussel numbers and quality as a function of the amount of fertilizer applied by Danish farmers (see Supplementary Methods) and winter and spring temperatures along the Wadden Sea coastline, using the same global climate database as described above (see “North-west Siberian climate”).

## Hierarchical modelling

We used an integrated hierarchical model 71, 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 71. In the following system, the state-space model, equation (1), expressing log-scale long-tailed duck dynamics (1968–2014) is written briefly as:

$$\begin{array}{c}{\widehat{n}}_{\left(t\right)}=a+c{n}_{\left(t-1\right)}+{\beta }_{R}{\text{log}R}_{\left(t-1\right)}+\\ {\beta }_{D}\left[{\left({1-w}_{D}\right)DIN}_{\left(t-1\right)}+{w}_{D}\left({-1*DIP}_{\left(t-1\right)}\right)\right]+{\epsilon }_{\left(t\right)}\end{array}$$

1

in which \({n}_{\left(t\right)}\) is the population size estimate of long-tailed ducks in year \(t\) that has a negative binomial error structure 49,72 around the expected population size \({\widehat{n}}_{\left(t\right)}.\) \({R}_{(t-1)}={1+p}_{\left(t-1\right)}\) is a variable for recruitment based on estimated juvenile proportion (\({p}_{\left(t-1\right)}\)), from equation (2), in autumn and winter preceding spring migration. Terms \(a\) and \({\epsilon }_{\left(t\right)}\) are intrinsic growth rate (intercept) and random environmental disturbance terms, respectively. \(c\)is a density-dependent parameter whose value can vary from zero to one. Dissolved nitrogen \({DIN}_{\left(t-1\right)}\) and phosphorus \({DIP}_{\left(t-1\right)}\) describe annual nutrient pools in the southern Baltic Sea (outputs from equation (4a)). The nutrient series implemented into equation (1) were scaled to have a mean of zero and variance of one. \({w}_{D}\) weights nutrient effects and \({\beta }_{R}\)quantifies the effect of juvenile proportion in autumn and winter on subsequent spring migration counts; these parameters are beta distributed varying between zero and one (uninformative prior beta-distribution with \(\alpha =1\) and \(\beta =1\)). \({\beta }_{D}\) is normally distributed regression coefficient. The observation model for equation (1) is specified with Gaussian errors and *a priori* assumed environmental disturbances (wind) as: , where \(\gamma\) is a parameter for east-west aspect winds \({x}_{\left(t\right)}\), and \(\tau\) is the standard deviation of a random observation-error process (Supplementary Methods). The joint-likelihood 71 of the integrated hierarchical model combining equations (1), (2) and (3) is summarized with Supplementary Table S1.

We constructed a logit model to generalize the variations of juvenile proportions during the period 1967–2017, based on the number of juveniles per adult female shot by Danish hunters:

$$\begin{array}{c}{\text{l}\text{o}\text{g}\text{i}\text{t} p}_{\left(t\right)}={a}_{p}+{b}_{L}{z}_{\left(t-\text{1,1}\right)}+{{b}_{C}C}_{\left(t-3\right)}+\sum _{k=1}^{2}{{b}_{\left(S,p\right)k}S}_{\left(t,k\right)}+{r}_{p\left(t\right)}\end{array}$$

2

Here, \({p}_{\left(t\right)}\) is the expected proportion of juveniles (in autumn and winter) affecting subsequent spring migration numbers (population size \({n}_{\left(t+1\right)}\), equation (1)) of long-tailed ducks. \({z}_{\left(t-1,j=1\right)}\) is an estimate for (log) lemming population size in the previous year based on the longest series of lemmings (Kara Sea, \(j=1\)) from equation (3). \({C}_{\left(t-3\right)}\) is a variable for winter climate (NAOI Dec−Feb) three years before. \({{b}_{\left(S,p\right)k}S}_{\left(t,k\right)}\) terms quantify the effects of precipitation \(\left(k=1\right)\) and temperature \(\left(k=2\right)\) in North-western Siberia in May and June on juvenile proportions, and these were weighted, *e.g.* for precipitation \(\left(P\right),\) as \({S}_{t,1}={w}_{\text{J}\text{u}\text{n}\text{e},1}\left(1-{P}_{\text{M}\text{a}\text{y}}\right)+{w}_{\text{J}\text{u}\text{n}\text{e},1}{P}_{\text{J}\text{u}\text{n}\text{e}}\). \({a}_{p}\) and \({b}_{X}\) are estimated regression coefficients, and \({r}_{p\left(t\right)}\) controls for random effects (full description in Supplementary Methods).

Lemming population dynamics \({z}_{\left(t,j\right)}\) estimated using three separate population time-series \(\left(j\right)\) of lemming abundances in the Western Taimyr Peninsula (1965–2017, \(t=0\) refers to 1965) was estimated based on the following state-space model:

$$\begin{array}{c}{\widehat{z}}_{\left(t,1\right)}={a}_{z}+{c}_{z}{z}_{\left(t-\text{1,1}\right)}+\sum _{k=1}^{2}{b}_{\left(C\right)k}{C}_{\left(t,k,1\right)}+{{X}}_{1}a+{{X}}_{2}b+{b}_{t }\left(t-1\right)+{ϵ}_{z\left(t,1\right)}\\ {\widehat{z}}_{\left(t,j=2|j=3\right)}={\widehat{z}}_{\left(t,1\right)}+{\delta }_{j}+\sum _{k=1}^{2}{\theta }_{k,j}{C}_{\left(t,k,j\right)}+{ϵ}_{z\left(t,j\right)}\end{array}$$

3

in which the error in lemming population size estimates \({z}_{\left(t,j\right)}\) following a negative binomial around the expectation \({\widehat{z}}_{\left(t,1\right)}\) (detailed in Supplementary Methods). \({a}_{z}\) levels population \(j=1\), and additive scaling parameters for populations 2 and 3 are denoted by \({\delta }_{2}\) and \({\delta }_{3} ({\delta }_{1}=0)\); similarly, \({b}_{\left(C\right)k}\) adjusts climate (\({C}_{\left(t,k,j\right)})\) effects for population \(j=1\), and scaling (additive to population 1) for populations 2 and 3 are denoted by \({\theta }_{k,2}\) and \({\theta }_{k,3} ({\theta }_{k,1}=0)\). Subscript \(k\) indicates precipitation in summer \(\left(t,k=1\right)\) and the previous autumn \(\left(t-1,k=2\right)\) for each population \(j\). \({b}_{t }\) is a parameter for a temporal trend effect. \({ϵ}_{z\left(t,j\right)}\) controls for random environmental disturbances. Parameter vectors \({a}\) and \({b}\) measure the consistency in three-year cyclic dynamics designed to dummy variable (0, 1) matrices \({{X}}_{1}\) (1965–1994) and \({{X}}_{2}\), (1995–2017) in which ones enable regular effects by three-year intervals (details in Supplementary Methods).

A log scale state-space model for fertilizer \(\left({F}_{t}\right)\) effect on nutrient pools measured in the Danish Straits \(\left(j=1\right)\) and the Baltic Proper \(\left(j=2\right)\), returning state estimates \({D}_{\left(t,k\right)}\) for DIN \(\left(k=1\right)\)and DIP \(\left(k=2\right)\) for the southern Baltic Sea in 1970–2016:

$$\begin{array}{c}{D}_{\left(t,k\right)}={a}_{k}+{c}_{k}{D}_{\left(t-1,k\right)}+{{b}_{k}F}_{\left(t-lag\right)}+{\text{s}\left(F\right)}_{t-lag,k}+{e}_{\left(t,k\right)}\end{array}$$

4a

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 effect on the nutrient pools. Non-linear effects are specified with a spline function for fertilizers generating two additive smoothing variables 73: \({\text{s}\left(F\right)}_{t,k}\) (Supplementary Methods). \({e}_{\left(t,k\right)}\) is a normal distributed random effect term (equation (4a)). Observed nutrient amounts link to the states via an observation model, equation (4b), with intercepts \({a}_{j,k}\) and standard deviation \({s}_{D\left(j,k\right)}\). The \(lag\) parameter is binomially distributed and can have values of 0, 1, 2 or 3 (years) (Supplementary Methods).

To estimate how biomass of mussels, the main food items for long-tailed ducks, 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 well-mixed 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:

$${\mu }_{\left(t,j=1\right)}={a}_{\mu }+{c}_{\mu }{\mu }_{\left(t-\text{1,1}\right)}+{{b}_{F}F}_{\left(t-{lag}_{\mu }\right)}+{\text{s}\left(F\right)}_{t-{lag}_{\mu }}+f\left(T\right)+{e}_{\mu \left(t\right)}+{\delta }_{\mu \left(t,1\right)}$$

5

where \({a}_{\mu }\) is intercept and \({c}_{\mu }\) is a density dependent parameter. \({e}_{\mu \left(t\right)}\) controls for state-process errors due to environmental variance and \({\delta }_{\mu \left(t,1\right)}\) for demographic variance 74,75. A function \({\text{s}\left(F\right)}_{t-{lag}_{\mu }}\) generates smoother 73 as explained for equation (4a), and \(f\left(T\right)\) parameterizes the weighted effects of mean winter temperature (December–February) on mussel biomasses two and three years later. The \({lag}_{\mu }\) can have values 1, 2 or 3 (*cf*. equation (4a)). Biomass estimates \({\mu }_{\left(t,j=2\right)}\) for Schleswig-Holstein follows the parameterization given in equation (5), excepting an independent demographic error term \({\delta }_{\mu \left(t,2\right)}\) and additive scaling parameter for an observation model (Supplementary Methods).

A logit model for flesh proportion in mussels \({f}_{\left(t\right)}\) for the period 1998–2013 was expressed as:

$${\begin{array}{c}logit\end{array}f}_{\left(t\right)}={a}_{f}+{\sum }_{k=1}^{2}{b}_{T,k}{T}_{k\left(t\right)}+{b}_{\left(f\right)F}{F}_{\left(t-{lag}_{f}\right)}+{e}_{f\left(t\right)}$$

6

in which the flesh proportion is predicted with parameters \({b}_{T,k}\) for winter \({T}_{k=1\left(t\right)}\) and spring \({T}_{k=2\left(t\right)}\) temperatures, and \({b}_{\left(f\right)F}\) for lagged fertilizer use \({F}_{\left(t-{lag}_{f}\right)}\). \({lag}_{f}\) can have values of 0, 1 or 2 (years). Estimates \({e}_{f\left(t\right)}\) control for random disturbances on flesh proportions (Supplementary Methods).