Combined Effects of Hydrological Drought and Reduced Food Availability on the Decline of the Little Penguins in South Australia

Droughts in many regions of the world are increasing in frequency and severity which, coupled with effects from anthropogenic water extraction and diversion, are reducing river discharges. Yet to date, few studies have investigated the impacts of hydrological droughts (i.e., reduced river outflows to the ocean) on seabirds. Here, we examined the consequences of the “Millennium Drought” on the local decline of an iconic Australian seabird, the little penguin (Eudyptula minor). We analysed monthly and annual penguin numbers in relation to river outflow, rainfall, the characteristics of the coastal waters (sea surface temperatures and chlorophyll-a concentrations), and local abundance of key predators and prey species. We found a negative association between monthly penguin numbers and both sea surface temperatures and river outflow. Annual penguin numbers were positively associated with southern garfish numbers (our local indicator of food availability) but negatively associated with annual chlorophyll-a concentrations. Our findings emphasizing the need for further research into the effect of hydrological droughts on seabird populations and for improved river management that account for potential downstream impacts on the coastal environment receiving freshwater from rivers.  

For example, upwelling processes, and resulting cold surface waters, may act differently on various prey species and thus may be beneficial to some seabirds (by increasing their prey availability) but increase mortality in other seabird species due to reduced food productivity (see Newton et al., 2009) or thermoregulatory stress. Understanding the impacts that climatic and oceanographic factors can have on seabird population decline is thus crucial but challenging.
Rainfall patterns are expected to change across the world as a result of climate change (Milly et al., 2005;Power et al., 2013). Droughts, in particular, are becoming more frequent and extreme in many regions of the world (Dai, 2013), resulting in reduced water outflow to the ocean and decreased oceanic productivity across multiple trophic levels (Lamberth et al., 2009;Lee et al., 2012;Arnell and Gosling, 2013;Auricht et al., 2018). Human interventions, such as water extraction and diversion, can greatly exacerbate climatic effects of drought, and even lead to hydrological droughts that are independent of climatic effects (Wada et al., 2013). Increased nutrient loading from freshwater outflows into the ocean stimulates phytoplankton and zooplankton production (Morgan et al., 2005), often leading to large aggregations of larval fishes and their predators near river estuaries (Cyrus and Blaber, 1987;Cyrus and Blaber, 1992;Grimes and Kingsford, 1996;Gillanders and Kingsford, 2002). Therefore, knowledge of the link between freshwater outflow and coastal trophic communities is critical but this relationship is often poorly understood (Gillanders and Kingsford, 2002;Gillson, 2011).
The Murray River is the largest river in Australia and a major source of freshwater entering coastal waters in South Australia. The Murray River discharges into two shallow and large freshwater lakes (Lakes Alexandrina and Albert also known as the "Lower Lakes", regulated by a series of barrages), mixing in a shallow estuarine-lagoon (the Coorong) before reaching the mouth of the river (thereafter referred to as the "Murray Mouth") ( Figure 1). Freshwater outflow into the coastal environment from the Murray River has declined significantly over the past century due to intensive drought periods and increased extraction for human use (Maheshwari et al., 1995;Thom et al., 2020). Human demands have reduced the average annual outflow at the Murray Mouth by 61%, meaning that water now ceases to flow to the mouth 40% of the time compared to 5% before increased extraction for human use occurred (Csiro, 2008). In addition, between 2001 and 2010, Australia experienced one of the worst droughts recorded since European settlement, known as the "Millennium Drought". This combined pressure resulted in reduction of both seawater exchange and freshwater outflows (Mosley et al., 2012;Geddes et al., 2016), as well as decreasing coastal productivity in the area (Auricht et al., 2018).
Little penguins (Eudyptula minor) are colonial seabirds that become central-place foragers during breeding, with most of their prey being captured within less than 60 km of their colony when provisioning chicks (Collins et al., 1999;Bool et al., 2007;Hoskins et al., 2008). As a consequence, their breeding success is strongly correlated to food availability (Chiaradia and Nisbet, 2006) and large colonies are expected to be located near areas of high productivity, such as freshwater inputs or upwelling (see also Fortescue, 1999;Poupart et al., 2017). Tracking studies at several colonies in Victoria have shown that many little penguins have a tendency to forage near estuarine environments or freshwater outflows (Collins et al., 1999;Hoskins et al., 2008;Preston et al., 2008;Kowalczyk et al., 2015). Kowalczyk et al. (2015) in particular found that little penguins within Port Phillip Bay (Victoria) foraged closer to the Yarra River mouth during drought years, but increased their foraging range in years of high rainfallwhen outflows from the Yarra River increasedlikely to follow the dispersed nutrients (as nutrients play a crucial role in regulating primary production and hence the food web; Mortimer et al., 1999;Cruzado et al., 2002;Sommer et al., 2002) and thus maximize resource acquisition.
Little penguins are generalist feeders that rely mainly on clupeidae fish species, such as anchovy (Engraulis spp.) and pilchards (Sardinops sagax) (Klomp and Wooller, 1988;Cullen et al., 1992;Chiaradia et al., 2003;Chiaradia et al., 2012). In South Australia, anchovies are the most frequently consumed prey species and account for 31-67% of the little penguin estimated prey biomass (Bool et al., 2007;Wiebkin, 2012). Both anchovies (Wedderburn and Suitor, 2012) and little penguins (Colombelli-Neǵrel, unpublished data) have been recorded in the Coorong and Lower Lakes system of the Murray River, and the Murray Mouth is well within the foraging range of the little penguins breeding in Encounter Bay (20-60 km; average 41 ± 4 km; Figure 1; Bool et al., 2007). Therefore, it is likely that some dependence between the Encounter Bay penguins and productivity of the Murray River estuary exists. This hypothesis is further supported by the fact that penguin numbers in Encounter Bay have experienced major declines in the early 2000s (Colombelli-Neǵrel and Kleindorfer, 2014;Colombelli-Neǵrel, 2015a), which coincides with the Millennium Drought (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). While many factors have been suggested to explain these declines in Encounter Baysuch as predation on land or at sea, low reproductive success or juvenile survival, as well as parasite abundance (Bool et al., 2007;Wiebkin, 2011;Colombelli-Neǵrel, 2015b;Colombelli-Neǵrel, 2015a;)only one study to date investigated the role of these suggested factors. Specifically, Colombelli-Neǵrel (2015a) reported that juvenile survival, but not reproductive success, was the most critical variable affecting little penguin population trends on Granite Island. Yet the relationship between penguin numbers in Encounter Bay and Murray River outflows has never been explored, which is the aim of this study. Based on the findings from Auricht et al. (2018), we hypothesised a significant and positive relationship between little penguin numbers and river outflow, via primary productivity and thus food availability (see conceptual model in Figure 2).

Study Area
Penguin numbers (yearly and daily numbers; see below: penguin numbers) were collected on Granite Island (Encounter Bay) and environmental data (rainfall, sea surface temperatures, chlorophyll-a concentrations, food availability, and predator abundance; see below: predictor variables) were collected around the coastal ocean zone between Granite Island and the Murray River (Figure 1; see also  (Wiebkin, 2011;. Since 2012, after the end of the Millennium Drought when a period of higher river flows occurred, the population stabilised to approximately 20-30 individuals (Colombelli-Neǵrel, 2020). All the other colonies within Encounter Bay (Wright Island, West Island, Seal Island and Pullen Island) are now considered extinct (Colombelli-Neǵrel and Kleindorfer, 2014). As its name suggests, Granite Island has a granite rock coastline. Little penguins come ashore at dusk to breed and moult in naturally excavated burrows and artificial nests (i.e., concrete tunnels and arrangement of rocks) located all around the island. The Murray Mouth is located approx. 20 km by sea from Granite Island (Figure 1).

Penguin Numbers
We used two datasets comprising daily and annual penguin numbers: the daily data span was January 2003 to December 2018; the annual data span was January 1981 to December 2018 (specific periods for which data were obtained are specified in brackets throughout the methods). Penguin numbers were obtained via (1) nightly counts (daily data [2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013][2014][2015][2016][2017][2018] and (2) annual population census (annual data 1981-2018). The nightly penguin counts report on the total number of adult little penguins counted by the Penguin Tour Guides on the North Shore of Granite Island during the Little Penguin Tours (see Figure 1). Tours were conducted for visitors to spot little penguins as they return from the ocean at sunset and ran all year. Additional nightly counts were conducted by the Flinders University research team in weeks when no penguin tours were scheduled. These additional counts by Flinders University were conducted under the guidance of one of the Penguin Tour Guides and followed the same guidelines as the Penguin Tours; hence the data collected by the two teams were comparable. Counts were conducted within two hours after dark and lasted approximately 1 to 1.5 hours. Population census data prior to 2001 were collected by R. Brandle over five sections of the island and interpolated for the entire colony by D. Colombelli-Neǵrel. Maps of these surveys confirmed that the area monitored before 2001 corresponded to half of the area surveyed since 2001. Surveys before 2001 established that the density of little penguin across the island did not vary until 2006. Therefore, the numbers obtained by R. Brandle were doubled for the interpolation. From 2001 onwards, annual population censuses were carried out each year (except in 2004) by a large team of local volunteers and a penguin ecologist. Each year, the island was divided into separate smaller sections and each section was thoroughly searched for presence or absence of penguin nests. Once a nest was found, its status was noted as active or inactive. A nest was recorded as active if it contained eggs, chicks, or adults, or had clear evidence of penguin presence such as fresh droppings, a strong penguin smell or recent nest excavation. A nest was recorded as inactive if none of the above criteria were found or if it had cobwebs at the entrance indicating that no adult penguin was regularly entering or exiting the nest. Nests were also marked with talcum powder to avoid double counting by different teams of volunteers. Censuses were conducted during the peak of the breeding season to ensure comparability of the data between years, and total population numbers were estimated as two adult penguins per active nest found. Little penguins are asynchronous breeders that can breed any time between April/May and January/February of the following year in South Australia (Wiebkin, 2012;Colombelli-Neǵrel, 2015a;Johnson and Colombelli-Neǵrel, 2020). As a result, the beginning and duration of each breeding season (and hence the peak of the breeding season) is highly variable among years (Allen et al., 2011;Wiebkin, 2012;Colombelli-Neǵrel, 2018;Johnson and Colombelli-Neǵrel, 2020). During our study, the peak of breeding on Granite Island shifted from winter/spring (~July-September) to spring/ summer (~October-December) following 2012 (Colombelli-Neǵrel, 2018).
The reported decreases in penguin numbers were not due to migration to other colonies based on mark-recapture data collected on Granite and West Islands between 1998 and 2003, which showed that only 0.9% of the 1255 marked individuals moved between the two islands, with no individual breeding on another island than their island of origin (see (Colombelli-Neǵrel, 2015a). Additional mark-recapture data in the area showed that little penguin adult survival rate in South Australia was less than 10% during this period (see also (Colombelli-Neǵrel, 2015a;Colombelli-Neǵrel et al., 2020). In addition, all adjacent colonies became extinct during the monitoring period (Wiebkin, 2011;Colombelli-Neǵrel and Kleindorfer, 2014).

Predictor Variables
The dynamics of marine productivity are intricate as many variables can interact (Behrenfeld et al., 2006). Therefore, to understand the relationship between penguin numbers and the Murray River outflow, we included other environmental variables in addition to river flow in our analyses to determine which variable(s) had the most influence on penguin numbers (see conceptual model in Figure 2). We specifically focused on rainfall, sea surface temperatures, and chlorophyll-a concentrations as these have been shown to influence the dispersion and abundance of fish species (Congdon et al., 2007;Nicolas et al., 2014;Syed and Ahmed, 2015;Auricht et al., 2018), and as a result, seabird breeding success and survival (e.g., Pinaud and Weimerskirch, 2002;Monticelli et al., 2007;Le Bohec et al., 2008;Cullen et al., 2009;Sidhu et al., 2012;Mills et al., 2020). In addition, chlorophyll-a concentrations are recognised indices of phytoplankton abundance and biomass (primary productivity) and often serve as indicators for prey availability (Monticelli et al., 2007;Lo-Yat et al., 2011). Similarly, sea surface temperatures are good indices of climatic variability as seabirds are known to be sensitive to extreme temperatures because of thermal stress (e.g., Hochscheid et al., 2002;Oswald and Arnold, 2012;Meyer, 2014) and often correlate with chlorophyll-a concentrations (Radiarta and Saitoh, 2008;Kavak and Karadogan, 2012). Rainfall, sea surface temperatures, and chlorophyll-a concentrations have all been found to correlate with survival in little penguins in other studies (Sidhu et al., 2012;Agnew et al., 2015;Ganendran et al., 2016). We present our conceptual model of how all the variables interact in Figure 2. While chlorophyll-a concentrations have a near immediate response to nutrients from river outflow (Black et al., 2016), there may be a time lag of up to 12 months between primary productivity and actual prey availability for marine top predators, such as seabirds (Price et al., 2020). Therefore, we include in Figure 2, and in our analyses, both the direct and time lagged relationships.
Daily outflow estimations (in mega-litres, ML; 1981-2018) from the Murray Mouth (based on the sum of estimated flows from all the barrages) were obtained from the Department of Environment, Water and Natural Resources (South Australia) (see also Auricht et al., 2018). Daily rainfall data (mm; 1981-2018) were obtained from the Australian Bureau of Meteorology database on the following website http://www.bom.gov.au/climate/data/using Victor Harbor meteorological station (distance: 2.8 km from FIGURE 2 | Conceptual model detailing the relationships between the environmental variables included our study (rainfall, sea surface temperatures, and chlorophylla concentrations), food availability (prey abundance and location), little penguin survival, and their predators (long-nosed fur seals). Straight arrows represent direct relationships, while dashed arrows represent time lagged relationships. We also present in grey other variables (upstream water extraction and diversion, other climate variables, ocean mixing processes and nutrient concentration) that were not included in our models but affect food availability and environmental variables.
Granite Island). Sea surface temperatures (°C; 1995-2018) and chlorophyll-a concentrations (mg/m³; 2002-2018) were only available monthly. Data on sea surface temperatures within a 20 km radius of ocean surrounding Granite Island were sourced from the Integrated Marine Observing System (IMOS; see http://imos. org.au/home.html). Chlorophyll-a concentration estimations used in this study were derived as monthly composite satellite imagery from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS), a sensor on board the polar-orbiting satellite 'Aqua'. MODIS-Aqua Level-3 chlorophyll-a products (at approximately 4 km spatial resolution) are free to access on the NASA Ocean Colour website (https://oceancolor.gsfc.nasa.gov). Although daily chlorophyll-a products are available from this source, cloud cover can result in large numbers of missing pixels. Instead, monthly composite imagery was used, where each image is an average of all clear observations over each calendar month. The average chlorophyll-a concertation was calculated in a region stretching 20 km south-west from the Murray Mouth (see Auricht et al., 2018 for additional details on the methodology), which captures the ocean region surrounding Granite Island.
For the annual data analysis, the environmental data were averaged for the three months prior to the moulting period as this is a critical period for little penguins' survival and includes a significant part of the breeding period (Ganendran et al., 2016). The beginning and duration of each moulting period on Granite Island varied among years but generally started in December-January and finished in March-April (Colombelli-Neǵrel, 2018). For our analyses, we adjusted this three months period based on the respective timing of the moult for each of our study years. We also included (1) the annual number of long-nosed fur seals (Arctocephalus forsteri) as an indicator of predation risk (Wiebkin, 2011) and (2) the annual amount of southern garfish (Hyporhamphus melanochir) commercially caught in the southern Gulf St Vincent as an indicator of food availability. Garfishes have been recorded in the diet of little penguins in other studies (Klomp and Wooller, 1988;Montague and Cullen, 1988;Cullen et al., 1992) and are the second most frequent prey item found in the diet of little penguins in South Australia (20%; Bool et al., 2007). Local and relevant data on their most frequent prey item (sardines) were not available. Long-nosed fur seal and southern garfish data were only available as annual datasets. Data on the annual southern garfish commercial catch (tonnes; 2007-2017) were obtained from the Fisheries Research and Development Corporation (https://fish.gov. au/report/191-Southern-Garfish-2018). Annual long-nosed fur seal population influence was also considered (recorded on nonconsecutive years between 1991 and 2018), using data representing the highest number of seals observed during the penguin breeding season on West Island, located 5 km from Granite Island (see Figure 1). These observations were recorded by a local community member during repeated boat surveys, and by the local rangers as part of regular surveys conducted for the Department for Environment and Water.

Statistical Analyses
All analyses were conducted in R 4.1.1 (R Core Team, 2019) using the packages Astsa, dplyr, forecast, fpp2, fpp3, ggplot2, GGally, imputeTS, MASS, timeDate and trend. We used standard time series regression (Shumway and Stoffer, 2017) and multiple hypothesis testing (i.e., hypotheses being tested simultaneously) to understand the associations between the time series of the predictors (rainfall, sea surface temperatures, chlorophyll-a concentrations, food availability, and predator abundance) and the decline of the little penguins (our outcome of interest), using daily (converted to monthly; see below) and annual penguin numbers. The main aim of a regression analysis is to investigate whether a change in an outcome variable can be explained by a set of selected predictors, while controlling for multiple potential confounding factors (Bhaskaran et al., 2013). Time series regression is being applied to estimate the effects of the predictors over time based on their effects on the previous period of time. Poisson regression could not be used as the variables were correlated, especially the predictor variables. Standard assumptions for multiple linear regression apply for time series regression (Hyndman and Athanasopoulos, 2018): the relationships between the variables should be linear, the error terms should be independent and identically normally distributed with a zero mean and a constant variance. To find the best model, we followed three steps: model specification (also referred to as model identification; during this step, we established the trends of the data, any correlation and seasonality, and selected a set of candidate models through variable or model selection process), model fitting (also referred to as model estimation; during this step, we estimated the parameters of the models and tested for the significance of the parameter estimates), and model diagnostics (during this step, we checked the assumptions of normality (QQ plot), constant variance (residual plots) and independence (ACF of residuals; the Auto-Correlation Function, which measures whether earlier data in the series have some relation to later data). Furthermore, we estimated the goodness-of-fit by comparing adjusted R 2 and the Akaike Information Criterion (AIC). The best model is the one with the lowest AIC and the highest R 2 ).
Missing data in time series regressions introduce bias. As the data collection was irregular in our study (see Methods and Figure 3), instead of removing incomplete records, we created a modified dataset for both the annual and daily penguin data using interpolation method following Moritz et al. (2015) as the time series plots of the variables showed either trend or seasonality. For the daily data, missing datapoints were estimated using an interpolation method, except for 2004 and 2011, for which too many penguin datapoints were missing to apply this method. Missing datapoints for 2004 and 2011 were therefore estimated using data from 2005 and 2012 respectively, as patterns observed in 2004 and 2011 were likely to be similar to those observed in 2005 and 2012. For the annual data, missing datapoints for the penguin (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)2004) or predictor data (long-nosed fur seals: prior to 1991, 1995-1998, 2001-2003, 2005, 2007; garfishes: prior to 2007, sea surface temperatures: prior to 1995; and chlorophyll-a concentrations: prior to 2002) were estimated using either interpolation or the moving average method. To test for the robustness of our interpolations, we used the forecast and imputeTS R packages (Hyndman and Khandakar, 2008; Moritz   Moritz and Bartz-Beielstein, 2017). The ImputeTS package is designed specifically to address imputation for univariate time series. We use the interpolation method following Moritz et al. (2015) as the time series plots of the variables showed either trend or seasonality. We first analysed the distribution of the missing vales using the plots for each variable (penguin data, sea surface temperatures, rainfall, chlorophyll-a concentrations, and river outflow) using the function plotNA.distribution within the ImputeTS package. We then did the imputation for univariate time series, plotted the results, and compared the functions na.interpolation and na.kalman within ImputeTS R package. We decided to use the results from na.interpolation as they were similar to na.kalman but simpler and faster computationally. We then aggregated the daily data (penguin count, river outflow and rainfall) to monthly data using the average function. Determining trends prior to building time series models is a critical first step for time series analyses. Therefore, we confirmed declining trends in penguin numbers in the time series with non-parametric Mann-Kendall trend tests (Mann, 1945;Kendall, 1975), which are used to check whether data collected over time have constant increasing or decreasing trends (i.e., Mann-Kendall trend tests check for difference in trends between earlier and later data). Statistics obtained from the Mann-Kendall trend tests are determined by the sequences and the ranks of the time series rather than the original values. Both datasets (monthly and annual) were then explored using exploratory data analysis and correlation analysis; variables, such as penguin counts and river outflow, exhibiting nonconstant variance were log-transformed (i.e., to stabilize the variance). After model fitting, we performed model diagnostics on the residuals by checking the ACF of the residuals to ensure independence of the data, along with normality and constant variance ( Table 1). Seasonality existed for the monthly rainfall and chlorophyll-a concentrations data; and such seasonality has been considered in the modelling. Variable selection method, the forward selection based on Akaike Information Criterion, was implemented to choose the subset of the predictors that was "best" in a given sense (Sheather, 2009). While we acknowledge that variable selection may affect the properties of the estimators, as well as the standard inferential procedures such as tests and confidence intervals (i.e., the regression coefficients obtained after variable selection are often biased and the P-values obtained from F-and t -statistics are generally much smaller than their true values; Sheather, 2009), no model selection is immune to these difficulties (Sheather, 2009). In this study, we used three criteria to evaluate the subsets of predictor variables: the Akaike's Information Criterion, the Deviance and the Residual Deviance (Burnham and Anderson, 2002). The corrected AIC, AICc (a bias corrected version of the AIC) is used for the annual data as the sample size was small. As stated above, we conducted two analyses for the annual penguin numbers: one including the direct relationships between the variables and one including the time lagged relationships described in Figure 2.

RESULTS
There was a clear reduction in daily river outflows from the Murray Mouth to the coastal ocean during the Millennium Drought, with minor discharges occurring between 2003 and 2006 and zero discharge between 2007 and 2010 ( Figure 3A). High river outflows started to occur again following 2011-2012 ( Figure 3A), which coincided with a stabilisation in the decline of little penguin numbers on Granite Island (Figure 4).
Chlorophyll-a concentrations in the ocean around the Murray Mouth and Granite Island were elevated in years with high river outflows and low in years when outflows were reduced ( Figures 3B, C). A comparison of the chlorophyll-a concentrations in a year of extreme-drought period with no river outflow (low chlorophyll-a concentrations; 2006) vsa period of high river outflow (high chlorophyll-a concentrations; 2012) is presented in Figure 3C. Monthly sea surface temperature and daily rainfall temporal patterns are presented in Figures 3D, E respectively. There was a strong and distinct declining trend for both the annual and monthly little penguin numbers over time (monthly data z = -10.56, N = 192, P < 0.0001; annual data: z = -5.70, N = 38, P < 0.0001; Figures 4A, B). For the monthly penguin numbers, while the data showed a small and significant increasing trend after 2010 (z = 7.65, N = 95, P < 0.0001; Figure 4A), the overall declining trend (comparing before and after 2010) was stronger than the small increasing trend after 2010 (z = -8.56, N = 97, P < 0.0001; Figure 4A). The best model explaining the monthly penguin numbers included river outflow and sea surface temperatures ( Table 2A). The model explained 39% of the variation in the data. Monthly penguin numbers were predicted to decrease by a multiplied factor of 0.94 for every 1 unit (°C) increase of sea surface temperatures (b = 0.94, 95% CI: 0.90; 0.97; Figure 5A) and by a multiplied factor of 0.86 for every 1 unit (ML) increase of river outflow (b = 0.86, 95% CI: 0.84, 0.89; Figure 5B). Analysing the direct relationships between the variables (see Figure 2), the best model explaining the decline of the annual penguin numbers included all predictor variables (river outflow, chlorophyll-a concentrations, rainfall, numbers of long-nosed fur seals and southern garfish commercial catch; Table 2B). The model explained 74% of the variation in the data and the strongest predictors were the numbers of southern garfish and chlorophyll-a concentrations. Annual penguin numbers were predicted to increase in the model by a factor of 1.4 [= exp (0.33)] with increasing of 1 unit (tones) of southern garfish numbers (b = 0.33, 95% CI: 0.17, 0.50; Figure 6A) but predicted to decrease by a factor of 0.37 [= exp (-1.01)] with increasing of 1 unit (mg/m³) of chlorophyll-a concentrations (b = -1.01, 95% CI: -1.65, -3.68; Figure 6B). Annual penguin numbers were also predicted to increase with increasing river outflow (b < 0.001, 95% CI: 7.44e-07, 1.68e-05) and increasing rainfall (b = 0.01, 95% CI: 0.002, 0.03), but to decrease with increasing long-nosed   fur seal numbers (b = -0.02, 95% CI: -0.003, -0.01). Adding the time lagged relationships described in Figure 2 did not change our results as the best model explaining the data remained model 3 ( Table 1). While the assumption of independence was satisfied in model 4 ( Table 1; with the time lagged relationships), the assumptions of constant variance and normality were not satisfied. In addition, the adjusted R 2 for model 4 was lower than for model 3 and the AIC and residual deviance of model 4 were higher than model 3 ( Table 2).

DISCUSSION
Seabirds around the world can serve as sentinels of environmental changes and help us understand how they are affecting wildlife (e.g., Wanless et al., 2007;Divoky et al., 2015;Thomsen and Green, 2019). In this study, we showed that (1) monthly little penguin numbers, an iconic and locally breeding Australian seabird, were negatively associated with sea surface temperatures and river outflow, and that (2) annual penguin numbers were positively associated with southern garfish numbers (our local indicator of food availability) but negatively associated with chlorophyll-a concentrations. The results of this study suggest that decisions regarding river water management should consider, not only human and terrestrial environmental requirements, but also the long-term impacts that this may have for the coastal environment outside the river system (see also Auricht et al., 2018;Hallett et al., 2018;Thom et al., 2020).
Chlorophyll-a concentrations are considered a reliable measure of primary productivity in an area (and hence food availability) because they are a proxy for phytoplankton production, which forms the base of the food chain as prey for zooplankton, fish larvae and other heterotrophic organisms (Monticelli et al., 2007;Lo-Yat et al., 2011;Lander et al., 2013). A previous study showed that decreasing Murray River outflow during the Millennium Drought had reduced primary productivity in adjacent coastal waters, and postulated there was likely consequences for higher trophic levels (Auricht et al., 2018). Specifically, the authors showed that primary productivity collapsed in the near-shore region during the Millennium Drought, and that the primary productivity plume was up to 60 km from the coast in years of high river outflows (Auricht et al., 2018), which would encompass Granite Island (Figure 1 and 3B).
In our study, we found a negative relationship between chlorophyll-a concentrations and annual penguin numbers, which contrasts with a previous study on little penguins showing that annual survival increased with increasing chlorophyll-a concentrations (Agnew et al., 2015). As previously stated, chlorophyll-a concentrations have a near immediate response to nutrients from river outflow (Black et al., 2016), while there may be a time lag of up to 12 months between primary productivity and actual prey availability for marine top predators, such as seabirds (Price et al., 2020). Yet additional analyses including time lagged relationships (see Results as well as Tables 1, 2) did not better explain the variation in penguin numbers. Interestingly, there was also a negative relationship between annual southern garfish numbers and chlorophyll-a concentrations. It is possible that other environmental variables, not accounted for in our study, may have impacted food availability (and hence this negative relationship) as small pelagic fishes (such as those preyed upon by little penguins) are often highly dependent on plankton-based food webs, which are very sensitive to short-term environmental changes in salinity (see Kelly et al., 2016). Furthermore, increased chlorophyll-a concentrations would increase the production of fish larvae, thereby also increasing the presence of predators and fishing intensity (see below) in the area, which in turn could have long-term negative impacts on little penguins. Finally, this negative relationship may be due to the low capacity of recovery of the Granite Island little penguin population as annual penguin numbers remained low and continued to decline after 2010 when chlorophyll-a concentrations started to increase. Indeed, previous study showed that the decline of the little penguins on Granite Island was linked to low juvenile survival (less than 2%; Colombelli-Neǵrel, 2015a). Delayed sexual maturity (in little penguins, individuals do not breed until they are two or three years old; Priddel et al., 2008) and a decrease in juvenile survival can add significant lagged effects on population size and recovery, especially for small population where birth and death rates can be extremely random (e.g., Legrendre et al., 1999). This suggests that the observed low juvenile survival may have had a stronger negative effect on little penguin numbers than any of the potential benefits gained from the increased river outflow (and consequent chlorophyll-a concentrations) once the population reached critical levels. Further investigation into the food web and association between South Australian little penguin numbers and chlorophyll-a concentrations are needed to fully understand their causal and temporal relationships. Monthly little penguin numbers were negatively associated with sea surface temperatures and river outflow. This aligns with studies in South Australia (Johnson and Colombelli-Neǵrel, 2020) and Western Australia (Cannell et al., 2012) showing that little penguins had a lower breeding success when sea surface temperatures were higher in the months before breeding. However, studies in Victoria showed that foraging effort and first year survival in little penguins were positively associated with sea surface temperatures (Sidhu et al., 2012;Berlincourt and Arnould, 2015) and little penguins had a greater breeding success when sea surface temperatures were higher in the autumn preceding spring breeding (Cullen et al., 2009). The negative association with river outflow also contrasts with our annual results, as annual little penguin numbers were positively associated with river outflow. Monthly penguin numbers would have been affected by the overall annual decline of the population but also the probability that individuals be present on land at certain time of the year. Hence the contrasting results may be explained if greater numbers of penguins were fishing at sea for longer periods during times of high river outflow. Results similar to this were found by Kowalczyk et al. (2015) in Victoria where penguins increased their foraging range when the outflow from the Yarra River increased.
Seals are recognised predators of penguins worldwide (e.g., Hofmeyr and Bester, 1971;Du Toit et al., 2004;Johnson et al., 2006;Lee et al., 2019;Bester et al., 2020). In this study, we found that annual little penguin numbers decreased when long-nosed fur seal numbers increased, supporting the idea that seal predation may have influenced the decline of little penguins in South Australia (Bool et al., 2007;Wiebkin, 2011). Yet the relationship between penguin and long nosed fur seal numbers was close to 1:1, indicating that a decrease in penguin numbers due to seal predation would be very slow. While the available data may have underestimated the number of seals in the study area, previous studies have shown that the relative importance of little penguins in long-nosed fur seal diets varies significantly across regions (Bool et al., 2007;Baylis and Nichols, 2009;Reinhold, 2014), suggesting that little penguin availability is not the main factor driving predation. It is therefore possible that, in years of low river outflow when primary production was low and food resources may have been scarce, long-nosed fur seals switched to preying on little penguins (as seen in Antarctic fur seals (Arctocephalus gazella) and chinstrap penguins (Pygoscelis antarctica) when krill abundance was low; Daneri et al., 2008). An increase in long-nosed fur seal numbers may also have resulted in reduced food availability for little penguins as seals can compete with penguins for food resources (Croll and Tershy, 1998;Lowther et al., 2020) and long-nosed fur seals also consume southern garfishes in South Australia (Reinhold, 2014). These hypotheses should be further tested with additional tracking and diet studies.
Our model predicted an increase in annual penguin numbers with increasing southern garfish numbers, indicating that a decrease in food availability may have been responsible for the annual decline of the little penguin populations in the area (on Granite Island but also on Wright, West, Seal and Pullen islands; see Colombelli-Neǵrel and Kleindorfer, 2014). A previous study showed that juvenile survival was the most important driver of little penguin population trends on Granite Island (Colombelli-Neǵrel, 2015a). While no information is available for South Australian little penguins, studies in Victoria have suggested that mortality of juveniles in their first years occurred mostly at sea (Reilly and Cullen, 1982;Dann, 1991), likely due to a combined effect of parasites and starvation (Harrigan, 1992). These studies and our results suggest that low juvenile survival on Granite Island (less than 2%; Colombelli-Neǵrel, 2015a) may be driven by low food availability (see also Reilly and Cullen, 1982;Dann, 1988), which in turn would drive the observed populations declines. A study conducted in 2006 on Granite and West islands showed that garfishes were the second most frequent prey item consumed by little penguins (Bool et al., 2007). Yet, the foraging behaviour and diet of South Australian little penguins is little studied with only two studies to date (Bool et al., 2007;Wiebkin, 2012), and additional dietary studies would be beneficial, especially considering that southern garfishes are a commercially fished species in South Australia. Further research is therefore warranted to determine if commercial catch of southern garfishes may have negatively impacted little penguins.
In conclusion, we found a strong association between little penguin numbers and southern garfish numbers (our local indicator of food availability), which may explain the local decline observed in penguin numbers in Encounter Bay during the Millennium Drought. The fact that the Granite Island little penguin population still has not recovered in 2020 (Colombelli-Neǵrel, 2020), despite larger outflows in 2012-2013 and at the end of 2016, suggests that the population may have reached some critical reduction in the number of breeding birds during the drought period. While this study focused on a single colony, further decrease in little penguin numbers (as recorded in South Australia; Colombelli-Neǵrel and Kleindorfer, 2014; Johnson and Colombelli-Neǵrel, 2020, this study) may be observed in the future because of climate change and/or reduced rainfall/river outflow, considering their high dependence with estuarine environments and freshwater outflows (Collins et al., 1999;Hoskins et al., 2008;Preston et al., 2008;Kowalczyk et al., 2015). Much of the research to date on seabirds and climate change have focused on the impacts of increased ocean temperatures (e.g., Sandvik et al., 2005;Chambers et al., 2011;Descamps et al., 2017;Johnson and Colombelli-Neǵrel, 2020). Similar to other studies on Arctic systems showing the importance of glacial outflows on marine productivity as well as seabird foraging and populations (Greḿillet et al., 2015;Urbanski et al., 2017;Bertrand et al., 2021; see also Michel et al., 2015), our findings highlight that other factors need to be taken into consideration, such as the potential importance of river outflow for the health and resilience of the coastal ecosystem, and our results should be considered in future river management strategies. Given the wide distribution of seabirds, their key role in ecosystems, and the fact that droughts are becoming more and more frequent (Dai, 2013), future studies, both within Australia and elsewhere, are needed to identify which species may be affected by hydrological droughts to further enhance both seabird conservation and river management.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by Flinders University Ethics Committee (E388-E449).