Black-legged kittiwakes as messengers of Atlantification in the Arctic

Climate warming is rapidly altering marine ecosystems towards a more temperate state on the European side of the Arctic. However, this “Atlantification” has rarely been confirmed, as long-term datasets on Arctic marine organisms are scarce. We present a 19-year time series (1982–2016) of diet samples from black-legged kittiwakes as an indicator of the changes in a high Arctic marine ecosystem (Kongsfjorden, Svalbard). Our results highlight a shift from Arctic prey dominance until 2006 to a more mixed diet with high contribution of Atlantic fishes. Capelin, an Atlantic species, dominated the diet composition in 2007, marking a shift in the food web. The occurrence of polar cod, a key Arctic fish species, positively correlated with sea ice index, whereas Atlantic species demonstrated the opposite correlation indicating that the diet shift was likely connected with recent climate warming. Kittiwakes, which gather available fish and zooplankton near the sea surface to feed their chicks, can act as messengers of ecosystem change. Changes in their diet reveal that the Kongsfjord system has drifted in an Atlantic direction over the last decade.


Extended Data
Other 14.6 (6.7) 1. 8 3.9 (3.6) 1. 8 Table 2. Two-sample permutation test results to compare wet weight contribution of Arctic and Atlantic species (See Table 2 Table 3. Correlation statistics between diet indices and explanatory variables (See Fig. 6 for summary). "Index" column indicates the corresponding diet index, "a" the intercept for linear regression using non-detrended values, "b" the slope, "N" the number of annual averages, "P" the p value, "R 2 " the coefficient of determination and "r" the Pearson correlation. Bold p values are significant.  Table 4. Linear regressions used to estimate fish lengths from otolith size together with references used to estimate age-groups of forage fish (Extended Data Figure 1).
Regression column represents the linear regression used (FL = fork length, TL = total length, SL = standard length, OL = otolith length and OW = otolith width), with range column indicating the length range (mm) of fish used for the linear models given in Source column.
References for age-group length estimates are given in the Age source column.

Supporting information
Text S1. Change-point analysis

Methods
The change-point analysis was conducted using frequency of occurrence data for Arctic and Atlantic species, as well as for explanatory variables (sea ice index, temperature, population size, clutch size and breeding success). Since the focus was on identifying possible step changes in the diet data, explanatory variables were constrained to exclude changes before 1997. Change-points are reported as the last year of the previous regime.
One step change for each dataset was identified using the cpt.mean function from the changepoint package 14 using the "AMOC" (At Most One Change) method 15 . The p-value for change-points was estimated using the maximum log likelihood method and was inverted to acquire a p-value, which interpretation is analogous to classical null-hypothesis testing 14,15 .
. These models emulated a situation with a gradual change of a similar direction complemented by a step change. Finally, 5) two linear regressions with a free slope (i.e. with an interaction term) and a change-point . This model imitated a situation where the time series demonstrated two different gradually changing regimes throughout years. Models for each variable were then compared and best fitting one was chosen using Akaike's Information Criterion (AIC). Only models 1 and 2 were used for explanatory variables that did not demonstrate a step change according to the change-point analysis.

Results
The change-point analysis indicated a shift in diet data and population size in 2006 (Table   S1). Similar change-point was identified for sea ice index in 2005, whereas temperature, clutch size and breeding success data demonstrated no significant change-points. Model selection confirmed the step change in Arctic and Atlantic species, but indicated that a linear model described sea ice index time-series better than piecewise regressions (Table S2)  wider dataset can be seen as a benefit for such application. Consequently, diet data from incubating birds and chicks were not removed. In this SI, we show that keeping the incubating birds and chicks in the dataset did not influence the results, and further that the simplified approach of using FOs instead of raw data with mixed models did not affect our conclusions.

Methods
Effects of breeding stage (incubating / rearing), bird age (adult / chick), and breeding colony (Blomstrand / Krykkjefjellet / Not recorded) on composition of binary kittiwake diet data (Arctic species, polar cod, T. libellula, Atlantic species, capelin and T. inermis) were analyzed using general linear mixed effect models (GLMM), logistic functions and year as a random effect with the lme4 package 16 for R 17 . The GLMM models were formulated as follows: was the binary column of diet data for 1982-2016, x the factor column specifying breeding stage, bird age or colony, and y the data frame. Similar effects on total regurgitate mass were analyzed using linear mixed effect models (LMM) and logarithm-transformed mass-based data for 1997-2016. The LMM models were formulated as log~ + Significance of these models was tested by comparing the models to a null model with removed x in the model formulation using Chi-square tests 18 . The significance was confirmed by estimating confidence intervals for model parameters using likelihood profile estimation.
Effects of breeding stage and bird age on environmental correlations (see Fig. 6) was examined by, both, GLMMs on binary composition 1982-2016 data and LMMs on logarithm-transformed mass-based 1997-2016 data. The models were constructed by using year as a random effect intercept and both stage and age as random effect slopes for each year.
The resulting model for binary data was formulated as was the binary column of diet data (Arctic species, polar cod or Atlantic species) and z the numeric column of an environmental factor (sea ice index, temperature or population size).
Marginal and conditional R 2 values were calculated from the models using the MuMIn package 19 and used to access how well models described the underlying data 20 .

Results
Wet weight was significantly affected by breeding stage with total regurgitate mass being 41 % heavier during chick rearing compared to incubation (Table S3). Further, samples collected from adults were on average 46% heavier compared to chicks. These total regurgitate mass differences were, however, not critical to overall conclusions in this paper, as we examined mostly the composition of diet samples.
Breeding stage and bird age significantly affected diet composition of the kittiwakes, whereas breeding colony did not (Table S3). Atlantic species were, in general, more common during chick-rearing and polar cod during incubation. These significant differences were, however, caused by a few years. Removal of incubating birds would have led to removal of 1984 samples altogether reducing contribution of Arctic species (Table S4) Table 3, Table S6). Both, binary and mass-based data on Arctic species significantly correlated with sea ice index (Table S6). This correlation was caused by polar cod. Further, binary data for Atlantic species correlated negatively with sea ice index and positively with temperature. Finally, total regurgitate mass correlated positively with sea ice index. It should be noted that while mixed effect models use the entire raw dataset and may illustrate the patterns in a dataset more precisely, they are also subject to more pitfalls than simple linear regressions using average values due to their complexity and larger dataset 21 . Mixed effect models may be "a more correct" way of describing patterns in the dataset, but their interpretation is also more difficult than the interpretation of correlations. Therefore we believe that the correlations and linear regressions in the main part of the article appeal to a wider audience than mixed effect models. Consequently, given the similar results from both methods, we opted for the simpler method in the main text of the article.      Table S6. Mixed effect model results on effects of explanatory variables on binary ("Bin") and wet weight ("Mass") diet data. "β" gives the slope of mixed effect models, "Min" and "Max" the 95% confidence intervals, "Np" the number of model parameters, "Df" the residual degrees of freedom, "χ 2 " the Chi-square test value, "P" the corresponding p value, R 2 m marginal R 2 value and R 2 c conditional R 2 value.