Pronounced long-term trends in year-round diet composition of the European shag Phalacrocorax aristotelis

Populations of marine top predators are exhibiting pronounced demographic changes due to alterations in prey availability and quality. Changes in diet composition is a key potential mechanism whereby alterations in prey availability can affect predator demography. Studies of long-term trends in diet have focused on the breeding season. However, long-term changes in non-breeding season diet is an important knowledge gap, since this is generally the most critical period of the year for the demography of marine top predators. In this study, we analysed 495,239 otoliths from 5888 regurgitated pellets collected throughout the annual cycle over three decades (1985–2014) from European shags Phalacrocorax aristotelis on the Isle of May, Scotland (56°11′N, 02°33′W). We identified dramatic reductions in the frequency of lesser sandeel Ammodytes marinus occurrence over the study, which was more pronounced during the non-breeding period (96% in 1988 to 45% in 2014), than the breeding period (91–67%). The relative numerical abundance of sandeel per pellet also reduced markedly (100–13% of all otoliths), with similar trends apparent during breeding and non-breeding periods. In contrast, the frequencies of Gadidae, Cottidae, Pleuronectidae and Gobiidae all increased, resulting in a doubling in annual prey richness from 6 prey types per year in 1988 to 12 in 2014. Our study demonstrates that the declining importance of the previously most prominent prey and marked increase in diet diversity is apparent throughout the annual cycle, suggesting that substantial temporal changes in prey populations have occurred, which may have important implications for seabird population dynamics.


Introduction
Marine environments are changing rapidly across the globe due to a range of anthropogenic activities, including pollution, overfishing and climate change (Halpern 2009;Poloczanska et al. 2013). These effects have altered the abundance and distribution of lower trophic organisms such as plankton, with consequences for mid-trophic level fish which are the principal prey for a guild of marine top predators (Cury et al. 2000;Frederiksen et al. 2006). Many marine top predator populations are declining markedly in association with these changes in prey availability and quality (Paleczny et al. 2015;Sydeman et al. 2015). Altered diet composition is a key potential mechanism whereby changes in prey availability can affect marine top predators (Reid and Croxall 2001;Cury et al. 2011). Several studies have demonstrated long-term changes in marine top predator diet, in particular seabirds (Miller and Sydeman 2004;Gaston and Elliott 2014). However, these studies have mainly been undertaken during restricted periods of the annual cycle, because of logistical challenges of obtaining diet data throughout the year. In seabirds, diet studies are usually conducted during the breeding season, from samples delivered by adults to offspring (Barrett et al. 2007). However, the non-breeding period is critically important for the population dynamics of seabirds, since most mortality occurs at this time (Weimerskirch 2002;Frederiksen et al. 2008). Thus, a key question in understanding the link between changes in prey availability and seabird population dynamics is the extent to which there have been long-term changes in non-breeding season diet, and whether these differ from those during the breeding season.
Our understanding of seabird diet outside the breeding period is largely based on indirect methods such as stable isotopes and fatty acid analysis (Owen et al. 2013;Kowalczyk et al. 2014) or samples from shot/dead birds (Blake 1984;Harris et al. 2015). Such studies have produced valuable insights into non-breeding diet, demonstrating marked differences from the breeding season, owing to a combination of altered prey availability (Kowalczyk et al. 2015), energetic constraints (Markones et al. 2010), habitat association (Ainley et al. 1996) and, in migratory species, altered locations (Ronconi et al. 2010). However, there is very limited information on long-term changes in non-breeding diet. Green et al. (2015) examined differences in breeding and non-breeding season diet in Cape gannets Morus capensis over a 30-year period. However, due to sporadic sampling, their trends analysis was restricted to the breeding period only. To our knowledge, no published studies have quantified long-term trends in non-breeding season diet composition in seabirds, and compared these with trends in breeding season diet from the same population.
In this paper, we analysed three decades of year-round diet in the European shag Phalacrocorax aristotelis (hereafter shag) collected on the Isle of May, south-east Scotland. The shag is a coastally distributed seabird that spends a large proportion of the day and every night on land (Wanless and Harris 1997). Full-grown shags regularly regurgitate pellets containing prey remains, which can be collected at accessible roosts, offering a rare opportunity to quantify year-round diet (Barrett et al. 2007). Shags show a flexible foraging strategy such that diet varies substantially across the species range. Lesser sandeel Ammodytes marinus (hereafter sandeel) is the dominant prey in many populations (Harris and Wanless 1993;Velando and Freire 1999;Lilliendahl and Solmundsson 2006), but at others, Gadoids (Gadidae), in particular saithe Pollachius virens, are the principal prey (Swann et al. 2008;Lorentsen et al. 2018). Seasonal variation in diet composition has been recorded in some populations in response to changes in prey availability (Velando and Freire 1999;Lilliendahl and Solmundsson 2006). Previous studies of the Isle of May population demonstrated that, in the late 1980s and early 1990s, the diet of shags consisted mainly of sandeels, with limited evidence of seasonal differences in diet composition Wanless 1991, 1993). However, the North Sea has warmed substantially over the past three decades (Høyer and Karagali 2016), which has resulted in changes in the distribution, abundance and diversity of many fish populations, including sandeel (Perry et al. 2005;Van Deurs et al. 2009;ter Hofstede et al. 2010). A recent analysis of diet based on prey delivered to shag chicks on the Isle of May demonstrated a marked decline in the proportion of sandeel, from 0.99 (1985) to 0.51 (2014), over this period, along with a concurrent dietary diversification (Howells et al. 2017). The authors attributed this dietary change to climate-mediated alterations in the availability of sandeels and alternative prey. Similarly, a community-scale analysis of seabird breeding diet at this colony demonstrated a decline in the importance of sandeels over the past three decades (Wanless et al. 2018). As local sandeel populations are resident (Boulcott et al. 2007), it is probable that any effect of environmental change on abundance or quality of these populations will affect both breeding and non-breeding diet of shags which over-winter on the Isle of May. Thus, we might predict a decline in the importance of sandeel in the diet throughout the annual cycle. However, sandeel availability varies among seasons since they are present in the water column during the spring and summer, but are buried in the sand during the winter, apart from a brief period when they emerge to spawn (Wright and Bailey 1993). Furthermore, environmental conditions, habitat use and energetic costs also vary between seasons Michelot et al. 2017). Thus, any changes in overall prey abundance or availability during the study might have different effects on diet composition at different times of the year. However, the issue whether long-term changes in diet composition outside the breeding season has matched trends observed in diet during the breeding season (Howells et al. 2017) is untested. Therefore, our specific aims were to: (1) quantify year round diet composition of shags over three decades; and (2) test whether dietary trends differ between the nonbreeding and breeding period.

Quantifying diet
The study was conducted between 1985 and 2014 at a European shag Phalacrocorax aristotelis (hereafter shag) colony on the Isle of May National Nature Reserve, Firth of Forth, south-east Scotland (56°11′N, 02°33′W). Shags are present on the island throughout the year, with a resident proportion of the breeding population joined in winter by migrants from other locations (Grist et al. 2014), allowing for the collection of pellets throughout the year. Pellets were collected opportunistically (mean number of sample days year −1 ± SD 23 ± 14; range 3-49) at roosts and breeding colonies using forceps, placed into a plastic bag and frozen. The breeding status and age of individuals that produced pellets was unknown. However, as chicks do not produce pellets, all samples were from full-grown (i.e. fledged) birds (Russell et al. 1995).
Samples were submerged in a saturated solution of biological washing powder (Biotex©) and heated at 40-50 °C, until all soft tissue and mucus was digested. Residual hard parts (e.g. fish otoliths, vertebrae and mouth parts, cephalapod beaks, mollusc shells and crustacea exoskeletons) were then identified to the lowest possible taxon using keys in Härkönen (1986) and Watt et al. (1997), allowing the presence/absence of each prey type to be recorded in each pellet. Sandeels Ammodytes spp. (principally, lesser sandeels A. marinus; Harris and Wanless 1991), the most frequent prey type recorded, have previously been classified in dietary studies on the basis of age (Harris and Wanless 1991;Howells et al. 2017). However, differentiating between sandeel age classes is generally not possible from otoliths obtained from pellets due to the effect of digestive erosion on otolith structure. Therefore, for the purposes of this study, all sandeels were aggregated into a single prey category. The presence of sand was also noted, since it may arise from accidental ingestion when foraging in sandy habitats and therefore be an index of prey species that live in these habitats, notably sandeels (Winslade 1974;Holland et al. 2005). The number of otoliths of each prey type in each pellet was then counted. Each fish has two otoliths, but due to the large numbers that may be encountered in a pellet and the potential for otoliths within a pair to undergo differential digestion, it was not possible to accurately match otoliths from the same fish. Therefore, each otolith was treated as an individual sample within each pellet.
Pellet analysis has been used to quantify diet in a range of seabirds, including shags, cormorants, skuas and terns (reviewed in Barrett et al. 2007). In appropriate study systems, large sample sizes may be obtained in a non-intrusive way throughout the year. However, quantifying diet from pellets involves two well-established limitations that must be considered when interpreting the data. First, due to differential rates of erosion, small or soft prey may be completely absent or under-represented in pellets, with larger prey, or those with more resilient body parts, more commonly retained (Barrett et al. 2007). For example, Johnstone et al. (1990) showed that in captive shags the recovery of otoliths from Sprat Sprattus sprattus, sandeel and Cod Gadus morhua was 17%, 20% and 52%, respectively. Accordingly, the most robust diet metric used to quantify prey in pellets is frequency of occurrence, in which items are scored on the basis of presence or absence. This method does not capture prey types that are completely digested, but accounts for any differential in digestion rates among prey types that are recorded by giving equal weighting to prey types irrespective of abundance in the sample. We also considered a second diet measure that is typically quantified from pellets, the numerical abundance of different prey types. This measure is more informative, but must be interpreted with care because it is more sensitive to the effects of differential digestion rates (Barrett et al. 2007).
A second limitation of quantifying diet from pellets is that the exact date when the prey were ingested is not known. However, the vast majority of pellets were fresh when collected, and they do not persist on rocks at our study colony because they disintegrate in rain or are consumed by herring gulls Larus argentatus, so we consider that pellets will have been produced within ca. 2 weeks of the sampling date.

Dietary response variables
For each pellet, we recorded the presence or absence of diagnostic remains (e.g. fish otolith, vertebra, bone, mollusc shell, cephalopod beak) of each prey type. Frequency of occurrence was then calculated as the percentage of pellets in which the prey type was found in each period within each study year. We focused our analysis on frequency of occurrence of the top five most abundant fish prey: sandeel Ammodytes spp., Gadidae (Cod Fishes), Cottidae (Cottids), Pleuronectidae (Flatfish) and Gobiidae (Gobies). All other prey types occurred in ≤ 10% of pellets and could thus not be analysed robustly, but due to their low prevalence in the diet, we consider the omission of these prey unlikely to significantly affect our interpretation of changes in diet composition.
Numerical abundance is typically quantified as the proportion of otoliths of a given fish prey type relative to all otoliths in the pellet. However, where the diet is dominated by a small number of prey types, as in this study (Sandeel 88% and Gadidae 7% of all otoliths), analysis of relative proportions leads to problems of interpretation, since a change in one prey type cannot be readily distinguished from a reciprocal change in the other. We therefore modelled number of sandeel otoliths relative to all prey otoliths and number of Gadidae otoliths relative to all non-sandeel prey otoliths. All other individual prey types occurred too infrequently for their relative abundance to be analysed. However, their summed contribution was < 5% of all otoliths.
Diet diversity was quantified by calculating sample-level prey richness, which was the number of prey types recorded in each pellet. Due to the effects of digestion on prey items, it was not generally possible to identify all body parts to species level, but to a higher taxonomic level which varied with prey type (fish: family; Crustacea and Mollusca: subphylum; Polychaeta: class). As prey richness is a count, the aggregate, annual prey richness (pooling all pellets in each year) was systematically higher than the sample average (sample-level prey richness: median 5; range 0-9; annual prey richness: median 12; range 6-14). However, as annual prey richness is a measure of the total number of prey types exploited each year, we included it in our analysis.

Defining breeding and non-breeding periods
For the purpose of this study, a study year commenced at the onset of breeding in one calendar year and ended at the commencement of breeding in the subsequent calendar year. To determine the timing of onset of breeding in each study year we calculated the month in which the population median egg laying date occurred, estimated from weekly observations at long-term monitoring plots : median day of year range 101-181; Newell et al. 2015;updated). In shags, average incubation duration of a clutch of three eggs, the modal clutch size in this population, is 36 days (Potts et al. 1980), with fledging occurring at a mean of 53 days after hatching (range 48-58, n = 35; Potts et al. 1980). Therefore, we defined each breeding season as the month of median egg laying date plus the following 3 months. This 4 month period was longer than the breeding period of individual pairs (~ 3 months), but was designed to capture the spread of laying that occurs in each year (Daunt et al. 2007). We found that 97% of all observations of breeding activity (defined as observations of incubating eggs or brooding chicks; n = 29,075) at our long-term monitoring plots occurred in this 4 month time window, confirming that it was a robust representation of the breeding period. The non-breeding period commenced in the first month after the breeding period until the last month before the month of median laying date in the following year (range of months: breeding: April-September; non-breeding: August-May; Supplementary Material Table S1).

Statistical analysis
All statistical analyses were conducted using the R programming software (version 3.4.0, R Development Core Team 2016). To test for temporal trends and effects of period (breeding vs non-breeding) on sample-level presence, relative numerical abundance and prey richness, we fitted Generalised Linear Mixed Models (hereafter GLMMs), using the 'glmer' function in the 'lme4' package (Bates et al. 2015). Binomial models with a logit-link function were fitted for presence and relative numerical abundance, and Poisson models with a log-link function for sample-level prey richness. For each of the sample-level dietary components we fitted a global model containing fixed effects of year, period and a year by period interaction. This framework allowed us to test for temporal trends, the differences between periods, and differing temporal trends between breeding and nonbreeding periods in each of the dietary components. Within each model, we also included random effects for month, year and month nested within year, to account for residual temporal autocorrelation. To account for overdispersion, we also included an individual, sample-level random effect in models of sandeel otoliths relative to all prey and Gadidae relative to non-sandeel prey (Harrison 2015). We did not consider sample date as an explanatory variable, since this variable had no clear biological relevance, due to the variable time elapsed between pellet production and collection.
To identify trends in annual prey richness, where there was just a single value per year, we fitted a Poisson GLMM with a log-link function. We subtracted 6 (the minimum annual prey richness value over the study) from each value, so that the data are consistent with the distributional properties of the Poisson distribution. However, we present the results and plots on the original, unadjusted scale. This step was not necessary with the sample-level prey richness data, as the minimum value was zero, i.e. pellets where no species were identified. Visual inspection indicated that the annual prey richness may be exhibiting non-linear trends. To test this, a global model containing both a linear and quadratic numeric fixed effect of year was fitted, along with a categorical, annual level random effect of year to account for overdispersion (Harrison 2015). We weighted each annual prey richness value by the number of pellets per year and included a fixed (offset) effect of log(number of pellets year −1 ) to account for any systematic change in annual prey richness with annual sample size.
In order to compare models with different fixed effects but the same random structure we used maximum likelihood in all models (Zuur et al. 2009). In each analysis, the fixed effect of year was centred on zero (by subtracting mean year from each value) and rescaled (by dividing the centred value by the standard deviation of year). The inclusion of all years in the analysis led to difficulties with model convergence. Preliminary analyses confirmed that this was caused by the inclusion of years where samples were not collected in both the breeding and non-breeding periods, so these were excluded from the modelling process (707 samples in 7 years; 1985-1987, 1994, 1998-1999, 2008).
Model selection was performed on the four models (null, year, period, and year by period interaction) for each variable using Akaike's information criterion corrected for small sample sizes (AICc), where the best-supported model was considered to have the lowest AICc value compared to alternative models. Models within two AICc (∆AICc < 2) of the top model were deemed as having similar levels of support (Burnham and Anderson 2002), unless they contained an additional parameter, in which case they were considered uninformative (Arnold 2010). Analysis was conducted according to an established protocol (Zuur et al. 2010), with the 'MuMIn' (Bartoń 2016) package used to obtain model selection outputs (see Supplementary Material for full details of model selection). Due to the large number of models, we only report those within 10 AICc points of the best model in the main text.
For figures and tables, annual means were calculated by pooling all samples in each period within a year. For presence, each mean value was calculated as the frequency of occurrence i.e. the percentage of samples in which the prey class was present. For numerical abundance, each mean value was calculated as the proportion of all otoliths of a given prey type relative to all otoliths. To aid comparison with frequency of occurrence, we converted numerical abundance proportions into percentages. Study years commenced at the onset of breeding, so each spanned two calendar years. All study years were retained in figures of annual mean data , with model plots presented over the range of years included in the analysis .

Temporal and seasonal changes in pellet composition
The best-supported model for sandeel presence contained an effect of year, period and a year by period interaction (Table 3; full model selection table presented in   Supplementary Material Table S3). Overall, sandeel frequency of occurrence decreased markedly in both the breeding and non-breeding periods. However, the decline was more pronounced during the non-breeding period, from 96% in 1988 to 45% in 2014, compared to 91% to 67% during the breeding season (data values: Fig. 1a; predicted values from model: Fig. 2a). The best-supported model for both Gadidae and Cottidae presence contained an effect of year only (Table 3; Table S3). Gadidae frequency of occurrence increased from 22% in 1988 to 66% in 2014 (data values: Fig. 1b; predicted values from model: Fig. 2b), whereas Cottidae frequency of occurrence increased from 5% in 1988 to 45% in 2014 (data values: Fig. 1c; predicted values from model: Fig. 2c; Table 3; Table S3). Overall, there was an increase in Pleuronectidae presence over the study, driven predominantly by the breeding period, when frequency of occurrence increased from 7% (1988) to 23% (2014), with frequency during the non-breeding period remaining relatively constant at 15% in 1988 and 14% in 2014 (data values: Fig. 1d; predicted values from model: Fig. 2d; Table 3; Table S3). Gobiidae presence increased overall between 1988 and 2014, but there was a significant interaction between year and period such that presence was higher during the non-breeding period at the start of the study (breeding 2%; non-breeding 6%), while by the end of the study the frequency was the same in both periods (breeding 21%; non-breeding 21%; data values: Fig. 1e; predicted values from model: Fig. 2e; Table 3; Table S3). Presence of sand displayed a substantial decline over the study, with a significant year by period interaction such that frequency reduced from 44 to 19% during breeding and 92 to 16% in the non-breeding period (data values: Fig. 1f; predicted values from model: Fig. 2f; Table 3; Table S3). Sandeel numerical abundance relative to all otoliths decreased from 100% in 1988 to 13% in 2014, but there was no evidence of a difference between the breeding and nonbreeding periods (data values: Fig. 3a; predicted values from model: Fig. 4a; Table 4; full model selection table presented  in Supplementary Material Table S4). The decline was less marked at the start of the study, but accelerated from the early 2000s. Gadidae numerical abundance relative to all non-sandeel otoliths reduced overall, but was consistently higher during breeding (data values: Fig. 3b; predicted values from model: Fig. 4b; Table 4; Table S4). The magnitude of change was similar in the two seasons, from 68% (1988) to 48% (2014) in the breeding period, and from 54% (1988) to 34% (2014) in the non-breeding period.  Sample-level prey richness increased over the study, but with a more marked increase during breeding (from 1.16 prey types pellet −1 in 1988 to 3.36 in 2014) than non-breeding (1.67 prey types pellet −1 in 1988 to 2.69 in 2014; data values: Fig. 5a; predicted values from model:  Table 5; full model selection table presented in  Supplementary Material Table S5). Annual prey richness displayed a quadratic trend over the study, increasing from 6.27 prey types year −1 in 1988 to 12.31 in 2014, with a peak of 15.80 in 2007 (data values: Fig. 5b; predicted Periods are reported as non-breeding (NB) relative to breeding. Table shows model rank compared to other models, model structure, fixed effect estimates, standard errors, z ratios, number of parameters (k), difference in AICc between top model and selected model (∆AICc) and Akaike weight relative to other models (ω i ). Due to the large number of prey types and models, we only report those models within 10AICc points of the top model, which is shown in bold (for full model selection tables see Table S3) values from model: Fig. 6b; Table 5; Table S5). However, a model containing a linear effect of year received similar support, providing strong evidence for an increasing trend in annual prey richness.

Discussion
We identified dramatic changes in the diet composition of full-grown European shags Phalacrocorax aristotelis  Periods are reported as non-breeding (NB) relative to breeding. Table shows model rank compared to other models, model structure, fixed effect estimates, standard errors, z ratios, number of parameters (k), difference in AICc between top model and top model (∆AICc) and Akaike weight relative to other models (ω i ). Due to the large number of prey types and models, we only report those models within 10 AICc points of the top model, which is shown in bold (for full model selection tables see  (hereafter shag) on the Isle of May over the past three decades both during and outside the breeding season. The dominance of lesser sandeels Ammodytes marinus (hereafter sandeel) decreased, with the decline in sandeel occurrence more marked during the non-breeding period. In contrast, the frequency of Gadidae, Cottidae, Pleuronectidae and Gobiidae increased. Prey richness also increased over the course of the study, in particular during the breeding period. These marked changes highlight the importance of monitoring changes in diet composition throughout the annual cycle.

Dietary change
Our findings of an overall decline in the dietary contribution of sandeel throughout the annual cycle, support our general prediction that changes in the importance of sandeels over time would be similar in breeding and non-breeding diets, since local sandeel populations are resident (Boulcott et al. 2007). One explanation for this year-round reduction is climate-mediated alterations in the abundance, availability or profitability of sandeels associated with rising temperatures in the North Sea (Arnott and Ruxton 2002;Van Deurs et al. 2009). Similar dietary changes have been observed in other seabird populations in response to changes in prey availability (Miller and Sydeman 2004;Gaston and Elliott 2014;Green et al. 2015). Howells et al. (2017) also recorded a reduction in the length of sandeels fed to nestling shags at this colony over the past three decades, which, due to the negative, non-linear relationship between calorific content and sandeel size (Hislop et al. 1991;Wanless et al. 2005), may be linked to the decreasing prevalence in shag diet. However, due to substantial digestive erosion of sandeel otoliths in pellets (Johnstone et al. 1990), it was not possible to  Periods are reported as non-breeding (NB) relative to breeding. Table shows model rank compared to other models, model structure, fixed effect estimates, standard errors, z ratios, number of parameters (k), difference in AICc between top model and top model (∆AICc) and Akaike weight relative to other models (ω i ). Due to the large number of prey types and models, we only report those models within 10 AICc points of the top model, which is shown in bold (for full model selection tables see use otolith length-fish length relationships to infer changes in sandeel length in this study. With flexible foraging behaviours, as evidenced by the wide range of prey types exploited throughout their range, shags may be able to adjust their diet in response to availability and quality of alternative prey. Such flexibility may be a key mechanism underpinning the dietary trends observed in this study, such that sandeel may have become scarcer or lessened in profitability compared to alternative prey, which may themselves have become more abundant or profitable. Data suggest that the energy density of alternative prey is similar to sandeels (Spitz et al. 2010). However, in the absence of estimates of prey availability or capture rates, it is not possible to fully establish the causes underpinning these temporal patterns in diet composition. Industrial fisheries may also reduce the availability of sandeels, with knock-on effects on seabird diet composition. However, the sandeel fishery off eastern Scotland did not overlap spatially with the foraging distribution of this shag population (Bogdanova et al. 2014). Furthermore, the fishery was only operational between 1990 and 1999 . As such, we would have expected a stepped reduction in sandeel occurrence in the diet over this period, which was not what we found. Similarly, Wanless et al. (2018) did not record a reduction in sandeel occurrence in the diet of the seabird community breeding at the colony during the 1990s. We therefore consider it unlikely that top-down fishing pressure was driving the observed trends in sandeel dietary contribution. The steeper decline in sandeel frequency of occurrence during the non-breeding period may be linked to reduced foraging capacity at this time of the year, as a result of shortened day length, adverse weather and absence of sandeels in the water column, apart from a brief period during spawning (Wright and Bailey 1993;Frederiksen et al. 2008;Daunt et al. 2014). Accordingly, any changes in overall prey availability over the course of the study might have had a more pronounced effect on diet composition at this time of year than during the breeding season. However, no seasonal difference in the rate of change was apparent in sandeel numerical abundance. This disparity with sandeel occurrence may arise because numerical abundance is quantified as the proportion relative to other prey, which themselves may have shown seasonal differences in trends. However, we could not test this since we could not distinguish changes in sandeels from reciprocal changes in other prey. Whatever the mechanism, the lack of difference between breeding and non-breeding periods in the trend in numerical abundance of sandeels relative to other prey suggests that this species has shown similar declines throughout the year in terms of biomass consumed. The overall reduction in frequency of sand is in line with these conclusions. Sand ingestion likely reflects accidental ingestion when foraging for sandeels, since shags generally extract sandeels directly from within the sand sediment (Watanuki et al. 2008), whereas other prey species that live in these habitats, such as Pleuronectidae and Callionymidae, are more likely captured on the sea floor.
The increase in dietary frequency of Gadidae accords with recent evidence of a distributional shift into Scottish waters of some Gadiformes in recent years (Cormon et al. 2014), including saithe Pollachius virens, the principle prey of shags is some populations. Pleuronectidae frequency also increased in the diet over the last 30 years, so shags may have continued to forage in sandy areas through the course of the study, but increasingly targeted Pleuronectidae, and other prey associated with sandeel habitats, such as Callionymidae, rather than sandeels. Gobiidae also increased, but this prey class is predominantly associated with rocky areas, which accords with past work on this population demonstrating the use of multiple habitats (Watanuki et al. 2008). Gadidae otoliths relative to other non-sandeel prey reduced over the study, suggesting that other non-sandeel prey have increased more rapidly than Gadidae. However, there was strong evidence that Gadidae numerical abundance relative to other non-sandeel prey was consistently higher during breeding. This is in contrast to Lilliendahl and Solmundsson (2006) who observed a higher prevalence of Gadidae in Icelandic shag pellets during winter. One possible explanation is that many Gadidae species use inshore waters as nursery grounds, with immatures moving into shallow, coastal feeding areas in the Firth of Forth during summer (Bergstad et al. 1987;Heessen et al. 2015).
One consequence of these dietary changes is that both sample-level and annual prey richness increased over the study, with the latter peaking in 2007. Long-term dietary diversification has also been observed in other seabird species in response to changes in prey availability (Gaston and Elliott 2014). The parallel increase in diversity at the single pellet and whole year scale suggests that, on average, the population is now exhibiting an individual generalist/ population generalist structure of resource use (Bolnick et al. 2003). Seasonal patterns of sample-level prey richness changed over the study, such that the increase was more pronounced during breeding, in line with seasonal differences in the pattern of change among Pleuronectidae and Gobiidae frequency of occurrence. Climate-mediated changes in fish populations have been widely reported in the North Sea, including changes in the abundance and distribution of many species (Perry et al. 2005;Dulvy et al. 2008). Thus, the dietary trends observed in our study population may be indicative of reductions in the abundance and availability of sandeel, increases in non-sandeel prey or a combination of both. These changes may vary among seasons, but without independent data on any abundance of these prey types it is currently not possible to distinguish these alternatives.

Limitations
It is important to recognise the limitations of estimating year-round diet from pellets when interpreting our results. The most important limitation of pellet analysis is the potential for underrepresentation of soft-bodied or easily digestible prey (Barrett et al. 2007). For example, Pholidae and Callionymidae (the otoliths of which are poorly sampled by pellet analysis) can form a substantial proportion of chick diet in this population (Howells et al. 2017), but were recorded infrequently in pellets. One important consequence of this is patterns of long-term change over time might have been different had we been able to detect all prey types. In particular, the increase in diversity over the course of the study may be greater than we could demonstrate if more digestible prey than sandeels have become more common in the diet throughout the year, as indicated from our diet data obtained from regurgitates (Howells et al. 2017). A further limitation of our study is that we had to pool all sandeel age-classes. As a result, we could not examine temporal and seasonal patterns in the relative contribution of different age classes, in contrast to our recent analysis of diet from regurgitations (Howells et al. 2017). Another consideration is that due to substantial differences in detection rates with sandeel size (i.e. larger fish are better represented in pellets; Johnstone et al. 1990), some of the observed reduction in sandeel relative numerical abundance may have been exacerbated by changes in detectability, since average sandeel length declined over the course of the study (Howells et al. 2017). However, given the dramatic trends observed in this study and the comparatively small decrease in sandeel size observed in chick diet (from in 104.5 mm 1988 to 92.0 in 2014), we consider our observation of a decline in sandeel abundance to be robust to this limitation. Finally, uncertainty in the date of pellet production could also have affected our results, for example by assigning pellets to the wrong period. However, given the length of non-breeding and breeding periods (several months) compared with the maximum likely duration between pellet production and collection (ca. 2 weeks), and the fitting of month as a random term in our models, we do not consider that this error would have had a strong impact on our results.

Demographic and conservation implications
The year-round reduction in the importance of sandeels in shag diet and associated dietary diversification may have important demographic consequences. In shags, the majority of mortality occurs in winter (Aebischer 1986;Harris and Wanless 1996;Frederiksen et al. 2008), linked to foraging capacity in more challenging environmental conditions (Daunt et al. 2006Lewis et al. 2015). Such changes may also be important during pre-breeding, when diet composition can be a key determinant of subsequent reproductive success (Sorensen et al. 2009). Prey availability during the breeding season is also a key determinant of breeding success (Daunt et al. 2001;Frederiksen et al. 2007). Crucially, effects on fitness are likely to depend on the relative profitability of different prey types throughout the annual cycle (Hislop et al. 1991;Litzow et al. 2004). Due to the difference in habitat associations between prey types, the dietary change observed may also have important implications for shag foraging distributions (Bogdanova et al. 2014;Michelot et al. 2017). The increase in proportion of non-sandeels in the diet could alter interactions with anthropogenic activities, such as offshore renewable developments or recreation. Shags in this population are partial migrants, whereby a proportion of individuals remain resident throughout the year while the remainder migrate (Grist et al. 2014). Studies that estimate diet composition during the non-breeding period throughout the population range would deliver a more complete picture of the potential implications for population dynamics and conservation management.
In summary, we identified substantial alterations in diet composition of a population of shags throughout the annual cycle over a 30-year period. Our results accord with recent climate-mediated changes in the distribution and abundance of many ecologically and commercially important fish species in the North Sea, most notably sandeel. To our knowledge, this study is the first to have quantified long-term trends in seabird diet outside the breeding season. The similarities and differences observed in these long-term trends compared with those during the breeding season highlight the importance of considering the diet of seabirds throughout the annual cycle in assessments of long-term dietary change. That the decline in sandeel frequency and abundance is apparent both during and outside the breeding season suggests that substantial temporal changes in prey populations have occurred, and may have important implications for seabird population dynamics in the region.