Estimating the breeding population of an elephant seal colony from a single census

Our aim was to develop a method for estimating the annual number of female elephant seals pupping in a colony from a single count. This is difficult because breeding females are not synchronous so there is no time when the entire population is present. We applied models that describe arrival and departure behavior to account for those missed in any one count and calculated correction factors that yield total population from any single count throughout a season. At seven colonies in California for which we had multiple female counts per year, we found consistent timing in arrival and departure both within and between sites for as long as 50 years. This meant that the optimal correction factor, the date when the maximum number of females was onshore, was consistent. At Point Reyes, a female count on 27 or 28 Jan can be multiplied by 1.15 to yield the total female population; at Año Nuevo Island, the correction was 1.17 on 25-26 Jan; and at Año Nuevo Mainland, 1.13 on 28-30 Jan. Correction factors at Channel Island colonies and King Range were also 1.13. Across the colonies studied, the factor 1.15 multiplied by a female count between 26 and 30 Jan was close to optimal, and this provides a method for estimating the female population size at colonies not yet studied. Our method can produce population estimates with minimum expenditure of time and resources and will facilitate monitoring of the elephant seal population size over its entire range.


Introduction
Northern elephant seals (Mirounga angustirostris) are large marine predators that aggregate on land at predictable times to raise pups, allowing straightforward counts of breeding females [ [1,2]]. There is never a day, however, when all females are ashore because some depart for all have arrived. Models are thus required to estimate arrival and departure curves and the total number of females breeding on the colony [ [3][4][5]]. Our objective here is to use these models to generate correction factors that convert a single count near the peak of the breeding season into an estimate of the total number of females using the colony in that year. We develop correction factors for seven colonies, then evaluate the consistency of census timing and resulting correction factors across sites. Accurate population estimates from single counts will allow efficient and inexpensive censusing of many colonies, especially valuable at sites that are difficult to access. This will contribute to long-term monitoring of the elephant seal pup production at all its major colonies from Mexico to Canada [ [6]].

Field sites
We used censuses of breeding females at seven colonies (Fig 1)

Female counts
A single datum is a complete count of all females onshore at one colony on one day. At Año Nuevo, we usually counted from high points on dunes near the females [ [2]], but counts from airplanes or drones were used to support ground counts on a few days when numbers were at their highest. Aerial counts confirmed that ground counts had no consistent bias. Most counts at Point Reyes were done from high cliffs overlooking the animals, allowing excellent visibility, but some were done on beaches immediately adjacent to the females [ [7]]. Counts at King Range were done from close proximity to the colony or from a trail within 50-100 m [ [8]]. At the Channel Islands, all counts were made from photographs taken from an airplane and have been evaluated for accuracy elsewhere [ [9,10]].

Census model
We define the census curve as the count of females through time during the breeding season. It is bell-shaped, starting from zero females in mid-December, peaking near the end of January, then dropping again until early March (Fig 2). Rothery and McCann [[3]] and Condit et al. [ [5]] developed models to describe the census curve precisely. The models yield estimates of daily arrival and departure curves, a(t) and d(t) (the total number of animals arriving [departing] on one day t), and cumulative arrival and departure curves, A(t) = t 1 a(i) and D(t) = t 1 d(i) (number of females arrived [departed]) through day t). The estimated census curve on day t is (1) The model parameters are estimated by fitting C(t) as closely as possible to the observed daily counts (Fig 2). Further details are presented in Condit et al. [[5]].
Here we add one substantial improvement. In earlier papers, we analyzed all years independently, estimating a set of parameters for each [ [2,5]]. Here, we fit the entire ensemble of counts at one colony in a single hierarchical model, with year as a random effect [ [11]]. Taking the mean arrival date as an example, each year has an estimated date, and those dates across all years are fitted to a Gaussian distribution known as the hyper-distribution, including one hyper-mean arrival date and one hyper-standard-deviation (Table 1). Each of the model parameters was constrained by a Gaussian hyper-distribution except total population. The advantage of the hierarchical method arises in years having too few counts to fit the full arrival-departure model. In previous papers, we omitted those years from calculations [ [2,5]]. With a multi-year model, they can be included because the arrival and departure curves from other years provide support (Fig 2).
The multi-year hierarchical model was fitted to five sets of data: the island and mainland colonies at Año Nuevo, Point Reyes, King Range, and the three Channel Islands combined; data were insufficient to model the Channel Islands separately. We did not attempt to create a single grand model with both colony and year levels, since there were too few colonies to support a random effect. All parameters of the model (Table 1) for those five datasets were estimated with a Bayesian Monte Carlo procedure, producing a posterior distribution and thus 95% credible intervals for parameters and the statistics calculated from parameters (for details on parameter-fitting, see references [5,12,13]). To compare colonies and to evaluate consistency of the census timing across years, we used the mean arrival date each year,â i , as fitted by the model, and the date of the maximum female census calculated from the fitted census curve C(t) in each year. Figure 2. Daily female counts through time and fitted census curves. Black points are observed counts. The black, bell-shaped curve is the fitted census curve, C(t); the blue dotted curve is the estimated cumulative arrival curve, A(t); and the dotted red curve the cumulative departure curve, D(t). A(t) and D(t) reach an asymptote at the total population of females using the colony in year i, Ni. Year 2016 at Año Nuevo Island (ANI) and 2019 at Point Reyes (PR) illustrate cases where the model generated standard curves even with counts too few to constrain the shape. It works because other years with ample counts set arrival curves via hyper-parameters µa and µs (Table 1). Other abbreviations: ANML, Año Nuevo Mainland; SNI, San Nicolas Island; KR, King Range. Range of the vertical axis varies greatly, but the dates on the horizontal axis are identical; vertical lines are at 1 Jan, 1 Feb, 1 Mar. Table 1. Six parameters for the census model in a single year and their hyper-parameters across years. The six main parameters are subscripted since there is one for every year i. The population size had no hyper-parameters, meaning every year's population was independent, and the two parameters describing tenure had hyper-parameters fixed in advance as narrow priors (µ d = 31.06, σ d = 0.27, µv = 0.1212, σv = 0.0062), a requirement for fitting the model [ [5]]. The remaining parameters had Gaussian hyper-distributions estimated from the data. In a model covering 20 years of censuses, 126 parameters are needed: 120 main parameters (6 per year) plus 6 hyper-parameters.

Hyper-parameter Parameter
Annual

The correction factor
Consider site s in year i and define N is as the total number of breeding females using the site that season. C is (t) is the estimated census curve at site s in year i, the number of females on the colony on each day t (Eq 1, Fig 2). Then C is (t)/N is is the proportion of females on the colony on day t. The inverse of the proportion is the correction factor: a multiplier for a female count on day t that yields the total female population, The correction factor is calculated separately in every year at a site, but to be useful in future years, we need the average across all years for site s, m is (t). To generate m is (t), census curves for each site s were simulated using many parameter combinations from the posterior distributions of the model parameters at that site; each simulation produced one m s (t), a multiplier for each day t at site s. To capture the total error associated with a multiplier, we incorporated both year-to-year variation and within-year error. For the former, we sampled 30 years at random from each site, and for the latter, made 40 random draws from the parameters' posterior distributions within each year; the mean and 95% quantiles of those 1200 simulations yield m is (t).
To be most relevant for the present, years were chosen from 2002-2018 at Año Nuevo Island, Año Nuevo Mainland, and Point Reyes. This was a period in which censuses and the model fit were consistent and reliable. There were only four years available at the Channel Islands and three at King Range, so all were utilized.

Census timing
Fitted censuses and observed counts followed a similar bell-shaped curve from late December to early March at all sites (Fig 2). Comparison of arrival date and peak census date across years and across sites demonstrates the consistency (Fig 3). In early years at Año Nuevo Mainland the arrival dates were around 17 Jan, but were consistently near 14 Jan after 1990. Point Reyes and Año Nuevo Island also advanced by 3-4 days in early years then stabilized ( Fig 3A). Between colonies, the long-term mean arrival date (the hyper-mean across years, µ a ) varied by 4.6 days: 10 Jan at the Channel Islands, 11 Jan at Año Nuevo Island, 14 Jan at Point Reyes, and 15 Jan at Año Nuevo Mainland and King Range (Fig 3A). The date of the peak census was 15-17 days after mean arrival: 25 Jan at the Channel Islands, 26 Jan at Año Nuevo Island, 28 Jan at Point Reyes, 30 Jan at King Range, and 31 Jan at Año Nuevo Mainland. Variation in the peak mirrored variation in arrival ( Fig 3B).
There were year-to-year fluctuations, however, mostly by < 2 days but with a few outliers (Fig 3). Some of the outliers were associated with unusual events. At Año Nuevo, 1983 had an enormous storm on 27 Jan, and surf covered the island beach. Most late arriving females used the mainland, creating early arrival on the island and late on the mainland [ [14]]. Again in 2010 at Año Nuevo and in 1995, 1998, and 2010 at Point Reyes, large storm surf forced females to move around, some to other colonies, thus creating unusual census patterns [ [15]]. Other outliers were associated with poor census coverage, for example in 1998-2000 at Año Nuevo [ [5]]. On the other hand, there was a steady delay in arrival over 2010-2015, consistent at four sites, supported by thorough census data and not associated with major storms. At the Channel Islands and King Range, the few years of data were consistent, but credible intervals were wide due to small samples (Figs 3). Correction factors for single counts Similarity in the shape and timing of the census curves at all sites meant similar correction factors (Fig 4). Using a female count at the peak of the census curve, the correction for yielding total population size, N i , was 1.13 at Año Nuevo Mainland, King Range, and the Channel Islands, 1.15 at Point Reyes, and 1.17 at Año Nuevo Island (Appendix Tables A1-A5). A count at the peak had the narrowest credible intervals. Within ±2 days of the peak, correction factors and credible intervals barely changed, but beyond ±5 days, they increased rapidly (Fig 4). The correction curve was early at Año Nuevo Island and the Channel Islands compared to the other three locations, corresponding to differences in female arrival and thus census timing. Nevertheless, at all five sites on 26-28 Jan, correction curves were close together, and the multiplier was between 1.13 and 1.19 (Fig 4).

Discussion
Over 50 years at Año Nuevo Island, the timing of the census curve has changed little. It was 4-5 days earlier than it is presently during early years, a pattern that repeated at Año Nuevo Mainland and Point Reyes, but after 1990 the census peak has been within 2 days of a long-term average in most years. The late census date in newly founded colonies is a consequence of dispersal behavior of immigrant females. At high-density colonies, the animals most likely to emigrate to new colonies are young females arriving late at a an existing colony, when density is high [ [16,17]]. This was conspicuous in 1983 at Año Nuevo, for example, when storms drove the latest females to the mainland [ [14]], and also at Point Reyes, which was established when the nearby Farallon Island colony was awash with storm surge [[18]]. In the closely related southern elephant seal, the peak of the female census curve is likewise consistent across colonies at the core of the latitudinal range [ [4]], but the outlying Peninsula Valdes colony, well to the north, was 10 days earlier [ [19]].
Consistency across sites produces consistent daily correction factors. At all sites we studied, the smallest error bars in corrections were between 23 and 30 Jan. Long-term consistency at Point Reyes and Año Nuevo suggests that the same week will continue to be optimal in the future. At the Channel Islands and King Range, where we do not have long time series, we cannot assess consistency, but for the present we assume it holds. Though colonies differed slightly in timing, the flat bottom of the correction curve leads us to suggest that 26-28 Jan is the best day for a single count at these locations and an optimal choice for counting other colonies not yet analyzed. Over those days, the correction factor was 1.13-1.19, and credible intervals mostly 1.1-1.2 (but as high as 1.3 at Año Nuevo Island). A count between 26 and 28 Jan thus produces a population estimate within 10% error, excellent precision for estimating population size of a highly pelagic marine predator.
The colonies we examined varied greatly in size. The Channel Islands each had over 10,000 females spread across many beaches, but at Año Nuevo, we separated the single beach on the island from several nearby mainland groups, and the King Range colony had only 100 animals. Despite these differences, the correction factors barely differed. We conclude that a population of breeding elephant seals can be estimated effectively from just one count, whether a single isolated group or a large number of groups counted together.
Correction factors are often applied to censuses of marine mammals because direct counts miss animals [ [20][21][22][23]]. In typical examples, animals come and go over a few hours, or sometimes days, and correction factors are strictly empirical. The case of elephant seals, with a very regular breeding haul-out, is the only one we know based on a quantitative description of behavior [ [3,5]].
A further approach at a novel location would be to use a few counts in the second half of January to establish the timing of the census and indicate which of the correction factors to use. Alternatively, if the observed peak count fell far outside 25-31 Jan, more work would be required to generate good estimators. It remains to be seen how closely the census timing holds at more distant colonies, such as those in Mexico [ [24]], which are five degrees latitude further south. In southern elephant seals, that is far enough for census timing to differ by 10 days [ [19]]. Once Mexican colonies are added to these analyses, the entire range of northern elephant seals -six sites in the United States and two in Mexico -would be covered [ [6,25]]. It would then be routinely possible, with single counts at eight colonies, to monitor the entire world population of a deep-ocean predator.