1.Introduction
South American sea lions Otaria flavescens (Shaw, 1800) are widely distributed, occurring along approximately 10,000 km of the South American Atlantic and Pacific coasts from Torres (Brazil) (^{Vaz-Ferreira, 1981}) to Cape Horn (Chile) and from Cape Horn to Zorritos (Peru) (Cappozzo, 2002). The species' distribution is characterised by features which allow breeding and foraging activities, such as solid substrate and food availability (^{Bowen et al., 2009}; Cappozzo and Perrin, 2009).
The pinnipeds' presence far from their usual breeding areas is often the result of erratic movements. However, some species like O. flavescens, display frequent and seasonal movements (^{Rosas et al., 1994}; Messias et al., 1994, Simões-Lopes, 1995). They are present in the Argentinean harbours of Mar del Plata (^{Rodríguez, 1990}) and Puerto Quequén (Westergaard et al., 1979) where they form stable colonies. In these haul-outs, individuals dwell most part of the year with displacement to rookeries only during the breeding season. Mating occurs from December to January with a few births in February. Hence, the breeding dynamics of the species is characterised by an intra-annual peak (^{Vaz-Ferreira, 1981}; ^{Campagna, 1985}).
Similar colonies are also found in two haul-outs on the coast of Rio Grande do Sul state (RS), Brazil: Ilha dos Lobos (29°20′ S, 52°06′ W) and Molhe Leste (32°11′ S, 52°04′ W) (^{Vaz-Ferreira, 1981}; ^{Rosas et al., 1994}) (Figure 1). It is presumed that individuals present in these sites come from the rookeries in Uruguay in search of food (^{Castello and Pinedo, 1977}).
Most sea lions found in both areas are adult and sub-adult males (^{Rosas et al., 1994}; Silva, 2002). Adult males, due to the absence of parental care, move away from the rookeries to return in the next breeding season (^{Rosas et al., 1994}). Females divide their time between foraging trips and fasting periods at nursing, therefore presenting more limited movements (^{Vaz-Ferreira, 1981}, ^{Campagna et al., 2001}). Preliminary results of ^{Giardino et al. (2009)} in a capture-recapture study suggested that males with high site-fidelity to the haul-out of Puerto Quequén in winter developed a particular strategy of mating that kept them connected with rookeries of Uruguay and Patagonia.
According to ^{Rosas et al. (1994)}, young individuals have less swimming skills and strength than older animals, which limit their movements to shorter distances and explain the small number of young animals (≤ 3 years of age) found in southern Brazil during winter months.
After heavy exploitation of sea lions in the first half of the twentieth century, protection measures aiming at population recovery were established. Meanwhile, the abundance of populations in Argentina is increasing (^{Dans et al., 2004}). However, population abundances in the Malvinas Islands and in Uruguay seem to be decreasing (^{Thompson et al., 2005}; Páez, 2005). Therefore, a better understanding on population dynamics is important to design effective conservation strategies.
Based on the possible Uruguayan origin of the sea lions occurring in southern Brazil and the availability of a considerable database of animal counting accumulated over a decade, we analysed the occupancy dynamics and abundance trends of O. flavescens in both haul-outs of Rio Grande do Sul, Brazil.
2.Materials and Methods
2.1.Study area and data collection
The coast of Rio Grande do Sul comprises a 618 km long stretch of sandy beach that is interjected only by the Tramandaí River mouth, on the north coast, and by the Patos Lagoon Estuary mouth, on the south coast (Figure 1).
The Patos Lagoon Estuary has a 20 km long entrance channel which is fixed by a 4 km long granite dike into the ocean (^{Seeliger et al., 1997}), known as “Molhe Leste (ML)”. The structure which was built to enable a safe navigation also became a suitable resting place for South American sea lions (^{Rosas et al., 1994}). The other site, named “Ilha dos Lobos (IL)”, is the only coastal island of Rio Grande do Sul deriving from the Mesozoic era basaltic spill. This island is located 2 km from the beach at Torres and has an area of 1,700 km^{2} comprising a levelled ground (^{Leinz and Amaral, 1966}, ^{Vieira, 1984}). Both areas are conservation units named “Refúgio da Vida Silvestre” (REVIS) and legally protected by Brazilian federal laws.
Between January 1993 and December 2002, 188 surveys were conducted comprising both locations. Monitoring activities consisted of one to four visits per month. Whenever possible, both sites were surveyed the same day and hour. Due to limitations in logistic and financial support, from ten years of monthly data collection, visits to ML and IL were not possible in 12 and 39 months, respectively.
In most pinniped populations, one cannot assume that all individuals are available at the time of local count (Buckland and York, 2002). If ignored, this fact will result in underestimation due to availability bias. However, since at ML and IL the maximum number of individuals occurs in the early morning hours, with a decrease during the day to increase again at dawn (^{Rosas et al., 1994}; ^{Sanfelice et al., 1996}), all field counts were set for 8:00 am.
Visits to both sites were made by boat, with the exception at ML on days with strong winds from the southern quadrant when the animals rested at the dike's central interspace, less visible from boats. On these occasions visits were made by walking the dike from land. We assumed same sighting efficiency for both approaches and used the data without distinction. Animals were counted by naked eye in one path transect with approximate duration of 15 min, simultaneously by two observers to ensure that animals were neither missed nor double counted. The animals in the water and near to the boat were counted as well.
2.2.Data analysis
We defined as the response variable y_{[i]} = “maximum number individuals observed at a given REVIS in month i ”, (i = 0, 1, 2,..., n), and considered it to follow a Poisson distribution with mean λ_{[i]}.
We further assumed a linear relation between log (λ_{[i]}) and the explanatory variables ‘Site’, ‘Month’, ‘Month^{2}’ and ‘Year’. The only categorical variable ‘Site’ was specified as:
x_{0[i]} =1, for “Molhe Leste” (ML)
x_{0[i]} = 0 for “Ilha dos Lobos” (IL)
The other covariates were numerical and, for computational convenience, rescaled to have mean zero and standard deviation equal to one.
The hierarchical structure of the model was as follows:
where β_{0} is the intercept and β_{1}, β_{2}, β_{n} are the coefficients of the explanatory variables.Since data y_{[i]} may show extra-variation characterised as an over-dispersed Poisson, we added a random effect with standard deviation σ to deal with this situation. When σ = 0 we fall back to a standard General Linear Model (GLM) or Poisson regression model. For σ > 0 we have a General Linear Mixed Models (GLMM) (^{Faraway, 2006}). Over-dispersion is often addressed by considering only fixed effects (σ = 0) but using a negative binomial distribution for y_{[i]} instead (^{Kinas et al., 2005}). However, the choice taken here is better suited to explore the intended hierarchical model structure (^{Royle and Dorazio, 2008}).
Twelve models with different combinations of explanatory variables, including interaction terms between ‘Month’ or ‘Month^{2}’ with ‘Site’, and with or without random effect were fitted to the data. Model parameters were estimated through Bayesian analysis (^{Gelman et al., 1995}) using Markov chain Monte Carlo (MCMC) algorithms to obtain the posterior distribution (^{McCarthy, 2007}). We used OpenBUGS in combination with R software through libraries R2WinBUGS and BRugs (^{Sturtz et al., 2005}; ^{Thomas et al., 2006}) to run the MCMC. Priors were chosen to be weakly informative. For each β_{j} this prior had a normal distribution with mean zero and variance equal to 10^{3}. The prior for the over-dispersion parameter σ was a uniform distribution between 0 and 10.
Models were compared by the Deviance Information Criterion (DIC) (^{Spiegelhalter et al., 2002}).
where , is the expected posterior deviance and p_{D} the effective number of estimated parameters. In hierarchical models the effective number of parameters can be much smaller than the actual number of parameters due to the build-in dependence structure among parameters. Furthermore, models for which differences in DIC and the DIC of the best fit (ΔDIC) are up to about 10 cannot be ranked with confidence (^{Spiegelhalter et al., 2002}). Therefore, to choose among models within this range of ΔDIC, considerations about the descriptive capacity of ecological important features should be taken into account as well. We used this approach here. A Bayesian residual analysis was also performed on the chosen model, based on the Pearson residuals (^{Dobson, 2002}). where y_{[i]} is the individual' counts and is the mean predict abundance.Within the Bayesian context we can obtain the posterior distribution for each Pearson residual and calculate the corresponding 95% posterior probability intervals. Only intervals which exclude zero indicate a bad fit. This objective residual check is unique to a Bayesian residual analysis and facilitates the identification of potential outliers.
Occupation probability ϕ_{[i]} is defined as the probability that at least one individual is present at month i (i.e. y_{[i]} > 0) and is therefore a function of the predict abundance λ_{[i]}:
Using this formula and the posterior distributions of each λ_{[i]} over the period 1993 to 2002, posterior distributions and 95% probability intervals for occupation probabilities were also obtained.
3.Results
The number of surveys per season is summarised in Table 1. The total number of counted individuals per survey varied from 0 to 121.
All posterior distributions were obtained with a Markov chain of 31,000 values, a burn-in period of 10,000 values and thinning of 2, resulting in a simulated sample of 10,500 values. Statistical diagnostic tools showed good convergence properties of the chains in all models and were omitted. Posterior means for the parameters in all proposed models are given in Table 2.
Models | Parameters | ||||||
---|---|---|---|---|---|---|---|
β_{0} | β_{1 }Month | β_{2 }Year | β_{3 }Month^{2} | σ Over-dispersion | DIC | ||
ML | IL | ||||||
1 | 3.282 | 3.242 | - | - | - | - | 5035 |
2 | 3.247 | 3.208 | - | 0.263 | - | - | 4701 |
3 | 3.211 | 3.166 | 0.386 | - | - | - | 4341 |
4 | 3.175 | 3.133 | 0.385 | 0.261 | - | - | 4015 |
5 | 3.546 | 3.526 | 0.534 | - | -0.435 | - | 3848 |
6 | 3.499 | 3.485 | 0.529 | 0.252 | -0.422 | - | 3552 |
7 | 3.241 | 3.841 | 0.371 (ML)1.010 (IL) | 0.255 | -0.061 (ML)-1.217 (IL) | - | 2936 |
1d | 2.777 | 2.523 | - | - | - | 1.334 | 1210 |
3d | 2.767 | 2.501 | 0.697 | - | - | 1.136 | 1209 |
4d | 2.762 | 2.519 | 0.694 | 0.424 | - | 1.059 | 1209 |
6d | 3.382 | 3.157 | 0.728 | 0.397 | -0.642 | 0.921 | 1208 |
7d(*) | 3.045 | 3.652 | 0.539(ML)1.012(IL) | 0.392 | -0.240(ML)-1.111(IL) | 0.802 | 1219 |
By adding covariates the DIC decreased indicating a better fit. However, strong downward changes in DIC values occurred consistently whenever, for a given model, the over-dispersion parameter σ was added, as can be seen when comparing models 1, 3, 4, 6 and 7 with 1d, 3d, 4d, 6d and 7d, respectively. This is a strong evidence that data are over-dispersed and that a mixed model is appropriate. According to DIC, the best fit (model 6d) included covariates ‘Month’, ‘Year’, ‘Month^{2}’ and over-dispersion σ. However, model 7d with month-site interaction includes important ecological features which are helpful to distinguish the abundance oscillation between sites. For being acceptably close to the best model (i.e: ΔDIC = 11) we chose model 7d for the remaining analysis.
The high posterior probability P(β_{0}[1] > β_{0}[2]) = 0.996 is strong evidence that the site-effect favours larger average numbers at IL in comparison to ML (Bayes factor (BF) = 249 indicating that this hypothesis is 249 times more probable than its alternative). In absolute values, this difference represents about 18 individuals (95% probability interval 5.5 to 32.6).
Since the quadratic term (Month^{2}) (Table 2) is smaller than zero (the posterior 95% probability intervals do not include zero) there is statistical evidence that abundance increases from January up to a maximum. This maximum is reached at different times in both sites, in August (winter) and October (early spring) for IL and ML, respectively. After reaching its intra-year maxima, abundance at both sites behave differently. While at IL it decreases, at ML high abundance is maintained until December.
The posterior 95% probability intervals (Pr.I.) of log(λ_{[i]}) obtained for model 7d and the corresponding observed log-counting log (λ_{[i]}) are further shown for ML and IL (Figure 2). The quadratic relation to abundance (Figure 3) and the marked seasonal trend (Figures 2), which is different for each site, are expressions of the Month-Site interaction.
A positive coefficient for ‘Year’ (Table 1, Figure 2) indicates a long-term increasing trend of abundance for both sites.
The residual analysis for model 7d is shown in Figure 4. Points indicate the residuals and vertical lines are the 95% Pr.I. Intervals above the line indicate months with model underestimates, while intervals below the line indicate model overestimates. All remaining points are within acceptable error range. Hence, for 86.1% of the data the model provides a good description of abundance, with overestimation in 11.1% and underestimation in the remaining 2.8% of the data.
Finally, the estimated occupancy probabilities for model 7d (Figure 5) define periods of the year in which area use was certain (probability 1). But, the increase over the years, in the smallest occupancy probability is a very clear indication that these sites are becoming more intensely occupied year after year.
4.Discussion
4.1.Covariates
The success of data analysis and inference depends of the choice of a statistical model that best approximates the data structure and the underlying ecological process. The principle of parsimony states that a good model should thrive for the smallest number of parameters able to represent the relevant features of data (^{Burnham and Anderson, 2002}). The choice among covariates to incorporate is, therefore, an important part in model selection and DIC a useful criterion to choose from competing candidates.
As with other penalised likelihood criteria (e.g., Akaike Information Criterium or Bayesian Information Criterium), we caution that DIC is not intended for identification of the “correct” model, but rather designed to merely compare a collection of alternative formulations (Carlin, 2000). In ecology there is often a great deal of ambiguity about the process involved in determining the true value of the response. If a number of possible model structures fit the observed data similarly, then there will be considerable uncertainty about which model is in fact the best one. If a plausible alternative model structure results in predictions that are different from those of the chosen best model, then there is some risk involved in bravely ignoring the alternative (Wintle et al., 2003).
In our analysis all selected covariates contributed to the decrease in DIC indicating that they help to explain the abundance variation. But, only the inclusion of the over-dispersion parameter was able to cause a decisive reduction in DIC. Similar preference for an over-dispersed model fit was achieved by ^{Kinas et al. (2005)} using a negative binomial distribution of the response variable in a conventional generalized linear model (GLM).
Despite the fact that DIC indicated 6d as the best fit, we opted for choosing the sub-optimal model 7d to make our inferences because this model describes important aspects of the distinctions between both sites. Although slightly over-parameterised, we believe that the ecological features that emphasise such distinction between the haul-outs is a fair price to pay.
The hypothesis formulated by Silva (2002) that sea lions, when leaving Uruguay after the reproductive season first arrive at ML, divide themselves between the two sites in a later time and return to Uruguay after a final “pit-stop” at ML, is corroborated by the distinct monthly occupancy patterns in both haul-outs.
The residual analysis was presented in the form of probability intervals. This is a unique feature of the Bayesian approach to residual analysis and a helpful tool to easily identify unacceptable discrepancies. Since in a Poisson regression variances are proportional to means, different residuals have to be analysed distinctively. This explains why some small residuals are considered unacceptable (short probability interval) in comparison to others that are much larger but associated to wide probability intervals.
4.2.Seasonality
Observed counts are useful as relative indexes of population size and indicators of abundance patterns (^{Márquez et al., 2006}). Due to a possible availability bias, the resulting abundance estimates should be treated as “minimal average abundance”. But, provided that the number of animals absent from the haul-out at time of counting represents an approximately stable percentage of total abundance, seasonal variation and long-term trends should be analysed as if they were bias-free.
The seasonal movements observed in our data may be partially attributed to changes in the availability of food resources near the rookeries (^{Rosas et al., 1994}; Stern, 2002), since food distribution in the ocean may be heavily influenced by physical properties like water temperature and nutrient availability (Bowen et al., 2002). The location of biologically productive sites and associated physical fronts may result in high abundance of prey that compound the diet of pinnipeds (Bowen et al., 2002).
The south Brazilian coast suffers influence of the Subtropical Convergence consisting of the confluence of tropical waters provided by the Brazil Current (BC) and sub-Antarctic waters provided by the Malvinas Current (MC). Wind-driven circulation can cause latitudinal displacement of this mixture zone (^{Palma et al., 2004}). These waters are also influenced by the continental water discharges of the La Plata and Patos Lagoon estuaries which support high biological productivity (^{Castro and Miranda, 1998}). All these factors are important features that impact fish abundance on the continental shelf, mainly demersal teleosts (family Sciaenidae) (^{Haimovici et al., 1997}), the preferred prey of O. flavescens.
Another important productivity source, ensuring high fish catches in the region, is the intrusion of nutrient through the South Atlantic Central Water along the bottom of the shelf reaching sub-superficial layers. This upwelling occurs mostly in the border region between the states Santa Catarina and Rio Grande do Sul, specifically at Torres where IL is located (^{Brandini, 1988}).
Furthermore, stochastic events like El Niño Southern Oscillation (ENSO) may carry out prey to the adjacent area of the Patos Lagoon estuary (^{Garcia et al., 2003}), causing the South American sea lions displacement in order to feed elsewhere. Thus, strong discharges caused by the ENSO may change the foraging behaviour including longer trip period and distances and increase of foraging effort (Bowen et al., 2002). This hypothesis corroborates the observations of Silva (2002) who described low abundance for O. flavescens in 1998 at ML, a year with strong El Niño. However, none of the covariates related to ENSO, like Southern Oscillation Index, Oceanic Niño Index and Multivariate ENSO Index had significant effect and were therefore omitted from the text.
Further influence on the foraging of pinnipeds is water temperature, which seems to condition the arrival of sea lions at the Rio Grande do Sul shelf. According to ^{Pinedo (1990)}, the displacement of O. flavescens is supported by the cold waters of the MC inner branch. The mean latitude at which the BC meets the MC depends on season (^{Reid et al., 1977}). The winter average latitude is 25.8° S and has been reported to reach the minimal latitude of 25.2° S in August (^{Souza and Robinson, 2004}). The Southern sea lions also reached the highest abundance at RS in winter and beginning of spring and specifically in August at IL which also agrees with the hypothesis put forward by ^{Pinedo (1990)}.
This study allows the conclusion that abundance of O. flavescens on haul-outs of Rio Grande do Sul varied seasonally reaching maximum numbers in winter and early spring with peaks in August and October for Molhe Leste and Ilha dos Lobos, respectively. This cyclic pattern is combined with a long-term linear increasing trend in abundance. This trend is further supported by the long-term increase observed for minimum intra-year occupation probability by O. flavescens. Finally, the hierarchical Bayesian inference of an over-dispersed Poisson model was able to provide a good fit to data and to contribute to better understand the interplay of various covariates in abundance and occupancy patterns.