There is more to climate than the North Atlantic Oscillation: a new perspective from climate dynamics to explain the variability in population growth rates of a long-lived seabird

)


Introduction
The past two decades have produced a plethora of ecological studies that hinge on the North Atlantic Oscillation (NAO) as a proxy for "climate" (Saetre et al., 1999;Ottersen et al., 2001;Forchhammer et al., 2002;Stenseth et al., 2003;Durant et al., 2004;Sandvik et al., 2005Sandvik et al., , 2012;;Dippner et al., 2014).This ubiquitous proxy usage may have led to the spread of a misconception of what NAO actually represents in an ecological context.Finding a statistical relationship without fully comprehending the underlying physics hinders the understanding of the climate-ecology interaction (Hallett et al., 2004;Stenseth and Mysterud, 2005).The point sometimes overlooked is that the NAO is just one of a number of modes of climate variability in the Northern Hemisphere (Barry and Carleton, 2001;Stenseth et al., 2003), as highlighted in Table 1.Thus, a better understanding of climate dynamics is critical before hinging on a proxy.
The NAO is defined as one of the leading modes of variability in the North Atlantic sector.It is characterized by an anomalous dipole of atmospheric pressure with centers over Iceland and Azores: a positive phase corresponds to a higher than normal subtropical high-pressure system and a deeper than normal lowpressure system called the Icelandic low (Hurrell, 1995;Barry and Carleton, 2001;Hurrell et al., 2003;Bader et al., 2011).The NAO could be considered as a manifestation of changes in the storm distribution: the cumulative effect of the storms and the associated changes in the storm track lead to the NAO signal (Mesquita et al., 2008(Mesquita et al., , 2011)).Nevertheless, the NAO is also influenced by other climatic processes, such as sea-ice changes in the Arctic and the Sea of Okhotsk in the North Pacific, and may shift its sign within a season (Bader et al., 2011;Mesquita et al., 2011).Even a single storm may have enough momentum to shift the NAO sign (Rivière and Orlanski, 2007) and the NAO is known to have a variable center of action (Zhang et al., 2008).
Indeed, with so many climatic features, it is no wonder the NAO "has recently become increasingly common to use" (Durant et al., 2004, p. 338).The use of the NAO as a "proxy, " without looking for other climatic clues, is implicitly based on the simplifying assumption that NAO is "climate variability."This may, in the worst case scenario, lead to the erroneous conclusion that a certain ecological process is unrelated to climate, only because no relationship with NAO has been found.Furthermore, Sandvik et al. (2012) have shown that the relationship between the NAO and population growth rate of seabirds is highly variable, both within and across species.This alone does not mean that populations that do not co-vary with the NAO are unaffected by climate, however.
Furthermore, when studies use NAO without providing a clear climate dynamics framework behind its effect on ecology, they become pointless (Mysterud et al., 2000).Statistical results need to be supported by physical mechanisms.Finding that there is a relationship with a proxy without providing a physical mechanism for it may even make the validity of the relationship dubious.Correspondingly, Ottersen et al. (2001, p. 9), in an early review of the ecological effects of the NAO, warn against an overuse of a proxy without understanding its full implications: "Unless a mechanistically based explanation is provided, the veracity of any statistical NAO-ecology relationship will thus remain uncertain."Moreover, relying on a specific proxy also makes it harder for ecologists to explain how atmospheric dynamics influence ecology (cf.Sydeman et al., 2014).
Thus, the objective of the present study is to illustrate standard statistical methods used in climatology to understand the connection between climate and ecological dynamics.Here, these methods will be applied to the population dynamics of the Common Guillemot Uria aalge, a long-lived seabird that breeds on Hornøya, a small island in northern Norway.Instead of starting with a proxy, this study will look one step before that, viz.finding what is affecting the species and where these climatic influences come from.Only afterwards, an appropriate climatic measure (or proxy) will be chosen.Our findings clearly demonstrate that there is more to climate than the NAO that impacts seabird populations.

Study Species, Study Area, and Data Collection
The Common Guillemot is a medium-sized (about 1000 g) seabird belonging to the auks (Alcidae).It is a colonial cliffbreeder with a circumpolar distribution in the boreo-low Arctic region.The Common Guillemot lays a single egg, which is incubated for 33 days on average.The chick is fed by both parents for 3 weeks, before it leaves the colony, still flightless, with its father.Birds start to breed when they are 5-7 years old, and annual adult survival is high (87-95%; Gaston and Jones, 1998).The Norwegian population amounts to approximately 15 000 breeding pairs, which represents a 90% reduction compared to the population size of the 1960s.The species is therefore listed as critically endangered in Norway (Kålås et al., 2010).
The fieldwork was carried out from 1980 to 2011 on Hornøya (70 • 23 ′ N, 31 • 9 ′ E), a 0.5 km 2 island in northeastern Norway.Annual counts of individual Common Guillemots on predefined monitoring plots were made by one of us (R.T.B.) late in the incubation period or during hatching.To minimize the day-today variation, five to ten counts were made on different days, and the mean number was used as an index of population size.The monitoring followed internationally standardized methods (Walsh et al., 1995).Successive annual estimates of the total population on Hornøya were based on a single, total count of 1900 individuals made in 1987 and the annual rates of change documented in the monitoring plots.See Barrett (1983Barrett ( , 2001) ) and Erikstad et al. (2013) for more details.
As can be seen from the population trajectory in Figure 1, there was a pronounced crash in population size from 1986to 1987(cf. Vader et al., 1990;;Erikstad et al., 2013).To estimate the yearly variation in population growth rate r, we used the variation in estimate of the total population size N in each census from 1 year to the next [r t = ln(N t /N t−1 )].

Mean Sea Level Pressure Data
Mean sea level pressure (hereafter "MSLP") from the European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis Interim Project (hereafter "ERA-Interim") was used throughout this study (Berrisford et al., 2009;Dee et al., 2011).

Extratropical teleconnection patterns Seasonal occurrence
Table adapted from Barry and Carleton (2001, pp. 396-397) showing the names, and acronyms when available, of the major modes of variability in the Northern Hemisphere together with their seasonal occurrence.Re-analysis data "provide a multivariate, spatially complete, and coherent record of the global atmospheric circulation, " and the ERA-Interim reanalysis, particularly, is one of the most sophisticated gridded datasets at high-resolution at present (Dee et al., 2011, p. 554).ERA-Interim data were created using the ECMWF Integrated Forecast System (IFS) model with fully coupled atmosphere, land surface, and ocean waves (Dee et al., 2011) together with the 4-dimensional variational assimilation (4D-Var) system (Courtier et al., 1994;Veerse and Thepaut, 1998;Dee et al., 2011).So, available observations from different sources (e.g.: satellite, station data) were assimilated into the model simulation, which allows for a product that is realistic (and in accordance with observational) data, with the advantage that it is on a three-dimensional lattice.In fact, rigorous quality control and experience from previous projects, such as ERA-40 (the precursor of ERA-Interim), have turned ERA-Interim into a state-of-the-art dataset for climate studies.ERA-Interim horizontal resolution is at T255 (nominally 0.703125 • , or around 79 km) and it contains 60 vertical model levels, with the highest being at 0.1 hPa.The horizontal resolution in ERA-Interim (79 km) is higher than that for another reanalysis product (with resolution on the order of 250 km) used in recent ecological studies (Dippner et al., 2014).
Thus, in the present study, 6-hourly ERA-Interim MSLP data from 1980 to 2011 were used.The winter season, referred to here as DJF (December, January, and February) was retained, starting from December 1979.These seasonal data were then aggregated into yearly averages per grid point.So, from here on, we will refer to the MSLP as starting from 1980(i.e.: December 1979, January and February 1980).This means, for instance, that when we discuss the year "1987" (outlier in the population growth rate data), we mean the winter 1986/87 (December 1986, January andFebruary 1987) for the MSLP variable.This provides the equivalent of 32 years of data, which conforms to a 30-year period considered by the World Meteorological Organizations to be the minimum requirement for climate studies (Parry et al., 2007;Baddour, 2011).

Wind Data for Polar Low Tracking
High-resolution data are needed to identify and track polar lows due to their small-scale compared with synoptic-scale lowpressure systems.For this end, another high-resolution data product is needed.Here, wind data at the 850-hPa atmospheric level from the NCEP Climate Forecast System Reanalysis (NCEP-CFSR) data were used (Saha et al., 2010(Saha et al., , 2013)).The NCEP-CFSR is a state-of-the-art reanalysis fully coupled product at the spectral T574 horizontal resolution (∼27 km) and 64 vertical levels, which is arguably one of the most advanced products at high-resolution to date.The importance of reanalysis data to climatologists was also expressed in Saha et al. (2010Saha et al. ( p. 1015)), "The general purpose of conducting reanalysis is to produce multiyear global state-of-the-art gridded representations of atmospheric states."Thus, the combination of an advanced storm-tracking algorithm (described below) with high-resolution advanced reanalysis dataset makes the result of this study unique.
Differently from ERA-Interim, the NCEP-CFSR was based on a fully coupled model system, the Climate Forecast System version 2 model (CFSv2), combined with an advanced data assimilation system (Saha et al., 2013); but like the ERA-Interim, NCEP-CFRS made use of the same observations from a number of sources in its assimilation system.These sources come from satellite data, aircraft and upper air balloon observations, as well as surface datasets.

Point Correlation Maps
Point correlation maps allow for the organization of large-data information and make it possible to identify climatic modes of variability (Wilks, 2011, pp. 67-70).This technique has been used in climate science for several decades.In fact, Wilks' (2011, p. 69) figure 3.28 shows the reproduction of a figure taken from the Norwegian climatologist Bjerknes (1969, p. 169).This figure illustrates that by using the method of point correlation for annual surface pressures around the globe with those in Djakarta, Indonesia, Bjerknes was able to identify a strong negative correlation with Easter Island, which reflects the atmospheric component of the El Niño-Southern Oscillation (ENSO).This is an example of a teleconnection pattern, that is, a remote phenomenon can affect other regions around the world (Glanz et al., 2009).It is interesting to note that Bjerknes' figure was also a reproduction of the earlier work of Berlage (1957), which exemplifies that this technique has been around for a long time as a means of understanding climate dynamics.
In this study, point correlation maps were created using the sample cross correlation at lag 0. Other lags could have been shown, but Sandvik et al. (2005 p. 823) find that the NAO has the highest correlation with the survival of Common Guillemot at lag 0 and here we test this effect.So, let X i,j be the vector of MSLP for latitude-longitude coordinate points i and j in the ERA-Interim grid, and Y be the vector of population growth rates for Common Guillemot.Thus, the sample correlation is determined as: where ρ i,j is the correlation coefficient at location i and j.Cov and Var represent the covariance and variance, respectively.The statistical significance of the correlation coefficient can be determined using Student's t-test (von Storch and Zwiers, 2002, pp. 77-78).

Point Regression Maps
Similarly, point regression maps can be created to aid the identification of hotspots where the predictor variable (e.g.: MSLP in this case) and the response variable (e.g.: population growth rate) are linearly related.These maps can be created through the least squares method, which is used to obtain meaningful information about the dependency relationships between variables (Draper and Smith, 1998).Also, the creation of such maps is often used in climatology and can provide further insights on how remote processes can affect other regions around a hemisphere or the globe (Bader et al., 2011;Mesquita et al., 2011).Regression maps are created by plotting the regression coefficient, that is, the slope (β 1 ) of the linear regression equation: where Y k;i,j , represents the k-th element of the response variable at location i and j, β 0 , and β 1 are the intercept and the slope parameters, respectively; and ε is the error variable.The slope is estimated β 1 as follows (Draper and Smith, 1998): where the overbar represents the mean values of X and Y at location i, j on the lattice.A t-test was used for significance for each grid point, which can then be overlayed on the map.A significance level of 0.1 was used to identify potential regions that can later be combined to construct an index.

Composite Analysis
Composite analysis is a well-known technique in climatology, which is used to construct climatic states that are typical (von Storch and Zwiers, 2002, p. 378).In order to illustrate the definition of composite analyses, we adapt the definition given in von Storch and Zwiers (2002, p. 378) as follows.Let z be a univariate index and M be the vector of the mean sea level pressure variable at a certain grid point.The composite is defined as: where represents sets of the index z and t represents time.The expectation of the composite is generally estimated as: for the observing times t 1 ,. . ., t k for which z t ∈ .In this study, composites were made for two sets: MSLP for the 5 years with the highest population growth rates ( h ) and MSLP for the 5 years with the lowest population growth rates ( l ).We estimate the difference between M h − M l applied to each grid point on the globe.In more general terms, composite analysis is the difference between selected states.For example, the difference between the mean of MSLP maps corresponding to years when the population growth rate of Common Guillemot is above/below a certain threshold with the mean of MSLP maps for years when the population growth rate is above/below a certain threshold average.This technique provides insights into atmospheric patterns that are more common for specific states of the population growth rate.

Storm Tracking Algorithm
In order to understand the 1987 crash, the storm-tracking algorithm TRACK, developed by Hodges (1995Hodges ( , 1999)), was used.TRACK is a computer algorithm, which has been applied in several storm-related studies (Mesquita et al., 2008(Mesquita et al., , 2011;;Bader et al., 2011;Zappa et al., 2014).Storms were identified in relative vorticity maps at the 850-hPa atmospheric level and were tracked temporally through a cost function.Relative vorticity represents the curl of the wind and how much it spins.It is determined from the u (west-east direction) and v (north-south direction) components of the wind, and is defined as (Holton, 2004, p.92): For the present study, polar lows, which are short-lived storms with a radius on the order of 100-500 km and surface wind speeds above 15 m s −1 , were tracked using the criteria detailed in Zappa et al. (2014Zappa et al. ( , p. 2603) ) with NCEP-CFSR (described above) as the input data.

Population Modeling
As a test of the utility of the different climatic indices in explaining Common Guillemot population dynamics, these indices were used as covariates in stochastic population models.Different models were compared using their AIC C (Akaike Information Criterion corrected for small sample size).All estimates are provided as means alongside their 95% confidence intervals.
Population dynamics of the Common Guillemot population on Hornøya were density-independent, as evidenced by the absence of a negative correlation between annual growth rates r t and population sizes N t (correlation coefficient R = +0.0005± 0.3664, p = 0.998), and the fact that the AIC C of the densitydependent (logistic) population model exceeded that of the density-independent (Brownian) model by 2.48 units.The model used thus had the form: where β i represents the slope of the ith environmental covariate X i ; ε, environmental noise, i.e., an independent variable with zero mean and variance σ 2 e ; N t , population size in year t; r, long-term intrinsic population growth rate; σ 2 d , demographic variance; X i,t , environmental covariate i in year t.The parameters β i , r, and σ 2 e were estimated using maximum likelihood (for details, see Saether et al., 2009;Sandvik et al., 2014), while σ 2 d was assumed to be 0.1, which is a realistic value for long-lived birds (Lande et al., 2003).To ensure that the latter assumption did not critically affect the results, different values of σ 2 d were tested; varying σ 2 d tenfold (0.01-1) changed the estimate of β by less than 2%.As the crash from 1986 to 1987 represented an influential outlier, population dynamics were analyzed with and without the population count for the year of 1986 included.

Correlation Maps
The point correlation between winter MSLP and the population growth rate variable are shown in Figure 2. Colors in red represent areas where the correlation is positive, meaning that high values of pressure in winter are associated with high values of the population growth rate that year (i.e., from the previous to next breeding season).An increase in pressure is related to a high-pressure system, which can promote cold air activity in winter, and a decrease in storm activity.Colors in blue are related to negative correlation between the variables.So in this case, a decrease in the values of winter pressure is associated with high values of the population growth rate.A decrease in pressure is related to low-pressure systems, which bring moisture and heat from lower to higher latitudes that create milder conditions in winter.
Figure 2A shows the correlation with the raw population data (includes the trend).It shows a clear dipole structure in the Arctic, as well as a wave-like propagation of positive-negative centers of correlation in mid-latitudes.The figure shows little or no resemblance to the NAO pattern in the Atlantic, but shows some resemblance, although shifted in center of action, to the Arctic Oscillation (Thompson and Wallace, 1998).It is, however, more similar to the Aleutian Low and Icelandic Low (AL-IL) seesaw pattern observed in later winter, another atmospheric mode of variability, which results in a dipole in the Arctic.The AL-IL is also the dominant pattern in late winter and it has consequences for climate in the surrounding regions (Honda et al., 2005).
The correlation with the detrended population data (Figure 2C) is quite similar in structure to that in Figure 2A.In fact, this pattern suggests that the trend in the population growth rate does not account for the seesaw-like pattern; it is a feature of another aspect of the population data.In order to investigate this, Figure 2B shows the data with trend, but removing the outlier in 1987 (i.e.: winter 1986/87 for MSLP; 1987 for r).Here, the features look very different indeed.The low pressure over where the colony is located is still there, but the seesaw pattern is much weaker.Also, for mid-latitudes, there is an increase in negative correlation.
Finally, the data with no trend together with the removal of the year 1987 is shown in Figure 2D.The pattern is similar to  Thus, 1987 is indeed an outlier, which seems to be related to a strong seesaw pattern in the Arctic.The main feature of the four figures is the low-pressure system over the Barents Sea and Pacific sector, as well as the high pressure over the Bering Sea.Both the outlier and the low-pressure system over the Barents Sea are investigated further in subsequent sections.But first, we investigate the aforementioned features through the use of point regression.

Regression Maps
Regression maps are given in Figure 3. Similar to correlation maps, the colors indicate regions of high (red) or low pressure (blue) associated with the population growth rate.Looking at the significant regions across Figures 3A-D, one can see some consistent patterns, such as the Bering Sea, Pacific Ocean and northwestern parts of North America, as well as the Barents Sea.If we focus on the latter, that is, the region where the colony is located, we observe that the regression with the detrended data without 1987 is significant.Thus, the trend and the outlier seem to mask the significance, showing that the pattern over the Barents Sea is the true underlying process (or signal) associated between pressure and population growth rate.
Regression maps are very useful in identifying "hotspots" of regions that have potential in explaining the association between the dependent and explanatory variables.Here they show, as pointed out before, that the pattern of the NAO does not emerge in Figures 2, 3.For instance, if a consistent dipole pattern between a region around Iceland and the Azores emerged, one would then relate it to the NAO.But this is not the case here.

Composite Maps
In Figure 4, composite plots of MSLP are displayed for the 5 years with the highest (i.e.: 1995, 1996, 1982, 2011, 2002) and the lowest (i.e.: 1987, 1985, 1984, 1988, 2007) detrended population growth rates.These plots provide further support to the robustness of the findings in Figures 2, 3, and are commonly used in climate studies.Figures 4A,B, show that the general patterns seem similar, however, the absolute values of the contour lines in the Pacific and Atlantic basin, as well as the distance between contours are different.The latter refers to the strength of the pressure gradient and provides an indication of the relative strength of the wind.So, in years with negative population growth rates, there is an intensification of the pressure gradient in both ocean basins compared to years with positive population growth rates.
The difference plot (Figure 4C) shows the difference between the MSLP for years with positive and negative population growth rates (i.e.: Figure A minus Figure B).Here, there is a striking similarity with the features observed in Figures 2, 3; more importantly, the dipole pattern in the Arctic between the Barents and Bering seas are evident.Thus, the fact that the feature observed in the Barents Sea is clearly seen in this composite confirms that the population growth rate is associated with the low-pressure system over the Barents region during winter.The dynamics behind this system will be explained further on.

Composite Analysis
Composite analysis has also been used to understand what the atmospheric features of 1987 were and how anomalous it was with respect to different climate baselines.Figures 4A,B, 5A based on the spacing of the isobars, Figure 5A would feature as the one with the least spacing between then.This is a first indication that the 1986/87 winter had an atmospheric feature that was stronger than that of the MSLP for the years with the lowest negative population growth rates.The bottom panel in Figure 5 shows anomaly plots, defined as the difference between MSLP for winter 1986/87 minus different winter baselines (from left to right): the mean of 1980 to 2011; the 5 years with the highest positive population growth rates in the detrended dataset (i.e.: 1985, 1984, 1988, 1998, 2007); and the 5 years with the lowest negative population growth rates in the detrended dataset excluding year 1987 (i.e.: 1982, 1992, 1995, 1996, 2011).These figures show a major correspondence with those in Figures 2-4, but with the colors shifted.This correspondence is especially seen for the Barents-Bering dipole pattern.The red color observed over the Barents Sea means that 1987 had an anomalous high-pressure system over the region, and as considered in the previous figures, a high-pressure system during winter in the Barents Sea is associated with a low population growth rate for that year.
Thus, the winter of 1986/1987 features as having an anomalous atmospheric configuration compared with different baselines.It had a major zone of high-pressure system, the socalled atmospheric blocking, over the Barents Sea.The latter feature was anomalous throughout the climatology of the region for the period considered here .

Effect of Local Weather
The results in the previous subsection indicate that the winter 1986/87 was indeed an outlier when it comes to the climatological conditions in the Barents Sea.It is therefore useful to look at the weather conditions that winter.Polar lows are a common feature in that region, which bring anomalously cold air from the pole, can have relative strength of a violent tropical storm (Nordeng and Rasmussen, 1992) and have damaging effects at the surface, such as icing (Icing implies the deposition and formation of ice on surfaces).Figure 6 shows mesocyclones, which are mediumsized storms, and polar lows in the North Atlantic and Barents Sea.The figure shows at least four polar lows crossing the central Barents Sea, which is the area where Common Guillemots from Hornøya spend the entire winter, according to ring recoveries (Barrett and Golovkin, 2000) and 3 years of geolocator data (Erikstad et al., unpubl. data).These polar lows were quite strong systems, as indicated by the vorticity, which represents how much the wind spins, with values higher than about 9 × 10 −5 s −1 .This suggests that through elevated thermoregulatory costs (due to low temperatures) and periodic difficulties in finding food (through high winds) the polar lows that winter could have increased the mortality of Common Guillemot.Indeed, many corpses of emaciated Common Guillemots were found washed ashore in Finnmark that winter (Vader et al., 1990).

Testing the Proxies
As a test of the utility of the indices described above, these indices, as well as the NAO index, were used as covariates in stochastic population models.When detrended, the two indices IDX2 and IDX3 improve the null population model (without covariates) by more than two AIC C units (Table 2).Both models are thus equally well supported, and considerably better so than the null model.Each of these two covariates explains 15% of the population dynamics.IDX2 and IDX3 with trends retained, as well as a pure trend model, are much poorer, in that they are roughly equally well supported as the null model (Table 2).NAO and IDX1 do not improve the null model, and explain less than 5% of the population dynamics.The finding that the NAO cannot physically explain the population dynamics of the Common Guillemot, is in accordance with the observation that NAO does not appear to be a mode of variability identified in Figures 2-5, either.
Models with two or more covariates did not improve upon the best models with one covariate (results not shown).Lagged covariates performed more poorly than their unlagged counterparts and/or the null model in all cases (time lags of 1-3 years were considered; results not shown).
The above results are based on population dynamics with the crash year removed.Even when 1986 is retained in the population time series, however, a detrended IDX3 is able to significantly improve upon the null population model (AIC C 2.5 units lower than null model, p = 0.026, R 2 = 0.148).The remaining indices (IDX1, IDX2, NAO) do not differ from the null model in this case (results not shown).

Discussion
Although the analysis of the population growth rate of Common Guillemot under classical climatological methods focused mainly on the mean sea level pressure variable, it addressed the main objective of this paper, that is, to highlight the need to understand the climate dynamics behind an ecological variable before constructing a so-called "proxy."In what follows, we present a discussion of the main results and propose a model to explain the atmospheric dynamical effect on the population growth rate variable.
The Hornøya population of the Common Guillemot, a species that is classified as "critically endangered" according to the Norwegian Red List (Kålås et al., 2010), declined by 80% in winter 1986/87 (Vader et al., 1990;Erikstad et al., 2013).The year 1987 was indeed an extreme event, as could be confirmed The models differ in the covariates included.Models are sorted in order of increasing AIC C (i.e., decreasing model fit).Models that are better supported than the null model (i.e., the model without any covariate) are highlighted using boldface.Estimates are provided with lower (LCI) and upper (UCI) 95% confidence intervals.The estimates of the trend model were multiplied by 10. by the anomalous high-pressure system over the Barents Sea, where the colony is located.It has been previously shown that the crash of Common Guillemot was related to low fish stocks (Vader et al., 1990;Erikstad et al., 2013), but it is possible that this situation was aggravated by (or even causally related to) the overlaying atmospheric condition in the winter 1986/87.From the clear and robust pressure system over the Barents region, as a response from the Pacific Ocean teleconnectivity, we could hypothesize that biotic responses (fish and/or birds) might have been forced by the atmospheric conditions over the Barents Sea.Note also that, the response in the Barents Sea is also forced through changes in sea ice, as shown in Figure 8.This highlights the importance of considering the climate system as a whole in the understanding of changes in ecological variables (cf.Sydeman et al., 2006).
There is mixed evidence for the relevance of NAO in ecology.NAO can explain variation in some different ecological variables in some species, but not in others (e.g.: Saetre et al., 1999;Forchhammer et al., 2002;Durant et al., 2004;Sandvik et al., 2005;Saether et al., 2009;Dippner et al., 2014;Sandvik et al., 2014).These studies and others illustrate the fact that NAO can be highly relevant for some variables and species; but an association may not mean causation.More is needed to assess NAO's usefulness, in this case, by finding a climate dynamics link to explain the association.The latter is clearly expressed in Forchhammer et al. (2002Forchhammer et al. ( , p. 1002) ) when discussing the impact of climate perturbations on arrival on breeding grounds and later survival and reproduction: "the timing of migration must be adapted ultimately to large-scale atmospheric systems."This large-scale referred to by the authors is the idea of looking at the bigger picture, before making a choice of an index to explain the variability.Therefore, looking for the dynamics first, before selecting an index, is here suggested as a more promising approach.Sandvik et al. (2005, p. 826) also recognize that: "causal pathways of these interactions are not easily identified.The NAO is defined as a temporal fluctuation of sea level pressure anomalies... and can as such hardly be said to be the cause of any of these responses.Rather it is convenient to treat the NAO as a "proxy" for different climatic processes...A full understanding of any climate-related response requires that one identifies the factors that cause it."The latter points to the need to look at the spatial contribution first to identify what affects the variability.In the present study, the variability explained by the NAO is not robust and does not appear in the spatial analysis of the relationship between MSLP and population growth rate, as concluded from the analysis of a long time series of the population growth rate of Common Guillemots.The local effects in the Barents Sea, which include the atmospheric forcing, mediated via teleconnectivity as found in this study, and the SST effects on the food web may act to affect the food web.This is also in line with the fact that there are positive correlations between SST and herring recruitment and stock biomass, which are relevant for Common Guillemot (Erikstad et al., 2013).
In the light of the results found in the present paper, and the discussions of previous scientific work, we propose the following mechanism to explain the variability in the population growth rate of Common Guillemot: • An anomalous low-pressure over the Barents Sea during winter, forced via teleconnectivity through the Pacific and local sea-ice changes, leads to the increase of moisture and heat brought by storm tracks into the region.This forces the SST locally and leads to warmer conditions-favorable conditions for Common Guillemot.• An anomalous high-pressure over the Barents Sea during winter, forced via teleconnectivity through the Pacific and local sea-ice changes, leads to blocking of storm tracks and loss of long-wave radiation from the surface to space.This leads to cold conditions both in the atmosphere and oceanunfavorable conditions for Common Guillemot.
Based on this model, the anomalous year of 1987 can also be explained.During the 1986/87 winter, the high-pressure was more anomalous than any other years in the study period.This led to extremely unfavorable conditions to Common Guillemot.With lack of food and increased energetic costs, these seabirds were at a great disadvantage.Any other external factor might be enough to cause a crash, such as the polar lows that brought even colder air from the pole and high winds into the region.We conclude that the role played by the Pacific as well as the Barents Sea is significant.It might indicate that the changes observed in the Barents Sea, where the colony is located, is triggered by a teleconnection from the Pacific via Alaska or via the north-east North America regions; these changes are also influenced by local changes in sea-ice.

Conclusion
The ubiquitous use of the NAO as a "proxy" in ecology and the search for alternative means of explaining variability have been analyzed in this research work.Although the NAO can explain the climate variability in the Atlantic sector, it is not the only teleconnection index available, and there needs to be a physical basis for the use of the NAO.So, the NAO as a "proxy" on its own should be used with caution.Also, before a decision is made on any proxies, one needs to see the bigger picture of what the possible climatic mechanisms are and how they can influence an ecological process.
Here, in search of such mechanisms, classical methods used in climatology to study climatic interactions were illustrated to study population growth rate of Common Guillemot: spatial maps of correlation and regression, compositing and a tracking algorithm to search for polar lows.We have used one climate covariate at the surface level for the analysis (i.e.: MSLP), and it has revealed an important climate mechanism at the Barents Sea.MSLP is just one of many other variables that could have been explored, such as: surface temperature, sea ice, temperature at different altitudes, geopotential height, precipitation, relative humidity, among others, to create a larger picture of the climate system and the interaction with the species studied here.Also, these variables could have been combined into a multivariate study and other advanced methods could have been used to explore their relationship.Thus, the present study illustrates the richness of information one gains through looking at ecology using the perspective that there is more to climate than the NAO.
Also, a possible physical mechanism was then proposed from the analysis, which linked teleconnectivity interaction via the Pacific into the Barents Sea, as well as interaction with sea ice changes in the Barents Sea region.An anomalous winter low-pressure system over the Barents Sea is associated with higher population growth rates, compared with years when an anomalous high pressure system is present: it creates warmer conditions brought by synoptic-scale storm track into the region.The opposite is true in anomalous winters with high-pressure system over the Barents Sea.The crash in 1987 could be explained by the extreme conditions of the severe anomalous high-pressure system over the region and the presence of polar lows.
One might ask if the techniques presented here apply to other species or variables.The answer is yes, but each species or variables would need to be analyzed on a case-to-case basis or through the use of multiple regression techniques.The latter is the subject of the next phase of this research work, where the idea of synchronization with different species and teleconnectivity will be explored.In doing so, other climate variables will also be addressed, such as sea surface temperature and geopotential height-in order to provide more insights into the dynamics of climate interactions and seabirds.
Overall, this research work has pointed to the importance of looking at what the data are telling in connection with climaterelated variables, instead of starting from on a specific "proxy."Only after the larger picture is obtained, can one narrow down existing modes of variability and/or create indices based on the response from the data.Thus, instead of looking at the data given the "proxy" (or mode of variability), [Data|Mode], one should explore the idea of looking at the climatic forcing given the data, [Forcing|Data], in search of modes of variability.Future work should address what is happening in the Pacific that triggers a wave-like pattern of teleconnectivity, as well as an integrated approach to understanding how climate affects seabirds, by looking at the role of climate forcing (ocean and atmosphere) on predator and prey in a multivariate fashion.
of Reading) for producing the tracks of the polar lows and Marine Ecology Progress series for the permission to reproduce Figure 1.We are very grateful to the National Seabird Monitoring Programme, which is funded by the Norwegian Environment Agency, for organizing the monitoring of Common Guillemots on Hornøya island.

FIGURE 1 |
FIGURE 1 | Dynamics of the Common Guillemot population on Hornøya, Northern Norway, from 1980 to 2011.Counts and estimates of numbers of breeding birds (filled circles, left-hand y-axis) and annual population growth rates r (open circles, right-hand y-axis).Figure adapted from Erikstad et al. (2013).

FIGURE 2 |
FIGURE 2 | Point correlation between mean sea level pressure (MSLP) and guillemot population growth rate (r) based on the period from 1980 to 2011.The four plots display the correlation coefficient

FIGURE 3 |
FIGURE 3 | Point linear regression coefficients between MSLP and guillemot population growth rate (r) based on the period 1980-2011.Each plot shows the regression coefficient between: (A) MSLP and raw r; (B) MSLP and r with 1987 removed from both

FIGURE 4 |
FIGURE 4 | MSLP Composite for: (A) Five years with the highest positive guillemot population growth rates (top left); (B) Five years with the highest negative population growth rates (bottom left); (C) Difference between 5 years with highest positive and 5 years with highest negative population growth rates.Units: hPa.
Figure 5 displays these composite maps.Figure 5A shows the MSLP configuration for the winter 1986/1987.Compared with Figure 4, the isobars (lines of same pressure level) seem to be much tighter both in the Pacific and the Atlantic basins.In fact, if one would rank

FIGURE 5 |
FIGURE 5 | Comparison between MSLP in 1987 with different baselines.(A) The MSLP in 1987.The bottom figures represent the MSLP in 1987 minus the MSLP in the following periods: (B) mean of 1980-2011;

FIGURE 7 |
FIGURE 7 | Regions selected for the MSLP-based index.

FIGURE 8 |
FIGURE 8 | Surface pressure regressed on sea-ice concentration in the Barents Sea.Sea-ice time series selected from the region indicated by the box.The data are based on the Coupled Model Intercomparison Project through the Max Planck Institute Earth System Model, CMIP5 MPI-ESM, control integration.Units: Coloring shows the change in hPa per one percent sea-ice concentration change multiplied by 10.

TABLE 1 |
Major teleconnection patterns in the Northern Hemisphere.

TABLE 2 |
Stochastic population models of the Common Guillemot population on Hornøya.